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2. 


SUMMARY 

An  important  impact  on  the  monitoring  of  a  CTBT  has  recently  been  made 
by  a  renewed  discussion  of  the  decoupling  problem.  Specifically,  since  the 
decoupling  of  an  underground  explosion  of  low  yield  by  detonation  in  a  cavity 
is  less  efficient  at  high  frequencies,  the  evasion  of  a  CTBT  through  decoupling 
becomes  increasingly  difficult  for  a  test  site  monitored  by  seismic  stations 
recording  frequencies  greater  than  10  Hz  over  high  Q  paths.  Of  particular 
importance  then  to  CTBT  monitoring  will  be  the  understanding  of  the  relative 
importance  of  scattering  versus  intrinsic  anelasticity  to  the  attenuation 
in  the  crust  and  lithosphere,  and  the  factors  that  cure  important  to  their 
regional  variation.  These  problems  are  treated  in  this  semi-annual  report  by 
comparing  the  predictions  of  multiple  scattering  theory  with  observations  of 
S  codas  recorded  from  earthquakes  in  the  Hindu  Kush  region  by  a  local  digital 


array .  The  complete  study  is  contained  in  the  Ph.D.  thesis  of  Ru-Shan  Wu, 
performed  under  the  supervision  of  Profs.  Aki  and  Toksttz .  The  work  included 
in  this  report  has  been  edited  from  that  thesis. 

The  study  has  examined  the  attenuation  of  local/regional  S  codas  from 
0.25  Hz  to  40  Hz  in  the  Hindu  Kush  region.  A  goal  of  the  study  was  to 
separate  the  relative  contribution  of  scattering  versus  intrinsic  anelasticity 
to  the  attenuation  of  coda  waves.  Coda  attenuation  has  been  analyzed  in  the 
frequency  domain  using  a  radiative  transfer  equation  technique,  which 
includes  multiple  scattering,  and  in  the  time  domain  using  weak  and  strong 
scattering  approximations.  The  frequency  band  less  than  1  Hz  appears  to  be 
dominated  by  strong,  multiple  scattering  and  wave  interference  that  cannot  be 
treated  using  the  radiative  transfer  equation  technique  in  the  frequency  domain. 

o''  S^TESTTflC  • 

air  !'r  r ;  ..  ■- 

NOT1. '  ,  ,-L2. 


TV.  '  ' 

o  ' 
Di 


J.  V  w  • 


^^•r.c,;:vnformntlonDivl.l°n 

Chief.  leot-.' 


v  .N  \\\'  \ 


This  band  may  be  dominated  by  scattering  of  surface  and  guided  mode  energy 
and  multiple  conversions  of  body  wave  energy  to  surface  wave  energy  and 
vice  versa  by  boundary  topography  and  elastic  heterogeneity.  In  the  band 
1.5  to  20  Hz  intrinsic  attenuation  dominates  scattering  attenuation.  The 
Hindu  Kush  data  indicate  that  this  intrinsic  attenuation  is  frequency 
dependent  in  the  1.5  to  20  Hz  band.  The  physical  mechanism  of  this  intrinsic 
attenuation  is  unknown.  In  the  band  above  20  Hz,  regional  phase  onsets 
cannot  be  easily  distinguished  and  scattering  attenuation  is  best 
described  by  diffusion  theory.  Although  coda  Q's  appear  to  be  similar  in 
tectonic  regions,  the  type  of  analysis  described  in  this  report  should  be 
applied  to  many  different  regions  before  general  conclusions  can  be  made 
about  the  relative  importance  of  scattering  and  intrinsic  attenuation  in 
different  frequency  bands. 
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Abstract 


In  order  to  separate  the  scattering  effect  from  the  intrinsic 
attenuation,  we  need  a  multiple  scattering  model  for  the  seismic  wave 
propagation  in  random  heterogeneous  media.  In  this  paper,  we  apply  the 
radiative  transfer  theory  to  seismic  wave  propagation  and  formulate  in 
frequency  domain  the  energy  density  distribution  in  space  for  a  point  source. 
Ve  consider  the  cases  of  isotropic  scattering  and  strong  forward  scattering. 
Some  numerical  examples  are  shown.  It  is  seen  that  the  energy  density  - 
distance  curves  have  quite  different  shapes  depending  on  the  values  of  medium 
seismic  albedo  B0  -  T|a/(ria+na) ,  where  qs  is  scattering  coefficient  and  T)a 
is  the  absorption  coefficient  of  the  medium.  For  high  albedo  (B>0.5)  medium, 
the  energy-distance  curve  is  of  arch  shape  and  the  position  of  the  peak  is  a 
function  of  extinction  coefficient  of  the  medium  qa  ■  tl*  +  rja.  Therefore  we 
can  separate  the  scattering  and  the  absorption  based  on  the  measured  energy 
density  distribution  curves. 

We  also  discuss  the  approximate  solutions  in  time  domain:  the  single 
scattering  approximation  and  the  diffusion  approximation.  We  apply  the 
formulas  of  diffusion  approximation  for  an  arbitrary  non-isotropic  scattering 
function  to  the  coda  envelope  and  discuss  its  relation  with  the  frequency 
domain  solution. 

The  data  from  the  digital  recordings  in  Hindu  Kush  region  are  used  as  an 
example  of  application  of  the  theory.  From  the  derived  energy  density 
distribution  curves  and  the  discussion  on  the  envelope  shapes  of  the  digitally 
filtered  seismograms,  we  conclude  that,  in  the  frequency  range  1.5  Hz  to  20 
Hs,  scattering  is  notthe  dominant  factor  in  the  measured  apparent 
attenuations,  i.e.  B0  <  0.5  in  the  Hindu  Kush  region  for  this  frequency  range. 
Due  to  the  insensitivity  of  the  shape  of  the  energy-distance  curve  for  the 
case  B0  <  0.5  end  the  fluctuations  of  the  data,  we  are  not  able  to  obtain  the 


Multiple  Scattering  and  Energy  Transfer  of  Seismic  Waves 
and  Application  of  the  Theory  to  Hindu  Kush  Region 
1.  INTRODUCTION 

Are  the  measured  apparent  attenuations  for  short  period  seismic  waves 
caused  by  anelasticity  of  the  media  or  by  scattering  of  the  heterogeneities  in 
the  media?  Is  the  single  backscattering  model  a  good  approximation  to  the 
coda  envelope  decay  or  do  we  need  a  multiple  scattering  model  which  will  have 
significant  differences  in  describing  the  coda  behavior  from  the  single 
backscattering  theory?  These  are  long-standing  problems.  In  order  to  answer 
these  questions,  we  need  to  develop  certain  multiple  scattering  model  for 
seismic  waves  and  compare  the  predictions  from  it  with  those  obtained  from  the 
single  scattering  theory.  O' Doherty  and  Anstey  (1971)  derived  a 
one-dimensional  multiple  scattering  formulae  for  a  stack  of  thin  layers  as 

|T(w)|  -  e-R(“)t  ,  (1.1) 

where  b>  is  the  angular  frequency  of  the  wave,  t  *  Nt  is  the  travel  time  of 
passing  through  the  stack,  r  is  the  travel  time  for  each  layer  and  N  is  the 
number  of  the  layers;  T(w)  is  the  transmission  response  and  R(u>)  is  the  power 
spectrum  of  the  reflection  coefficient  series  normalized  by  the  travel  time. 
The  exponential  form  of  (1.1)  itself  exhibits  the  indiscriminability  of  the 
multiple  scattering  effect  from  the  intrinsic  absorption,  if  we  observe  only 
the  decay  of  the  transmitted  waves.  Richards  and  Menke  (1983)  did  some 
numerical  experiments  on  this  model  and  discussed  some  possibilities  of  using 
the  relation  between  amplitude  spectra  and  phase  spectra,  the  frequency 
contents  of  the  coda  and  that  of  the  main  arrival  etc.  to  distinguish  the 
multiple  scattering  effects  of  thin  layers  from  the  intrinsic  attenuation.  We 
note  that  the  formulation  of  the  problem  by  O' Doherty  and  Anstey  is 
essentially  that  of  the  random  slab  problem  (see  Kay  and  Silverman  1958, 
Hoffman  1964).  The  results  are  presented  as  the  relations  of  transmitted  or 


reflected  waves  with  the  slab  thickness,  which  do  not  necessarily  represent 
the  aaplltude  attenuation  with  distance  or  the  envelope  decay  with  time  of 
seismic  waves. 

Kopnichev  (1977)  formulated  the  double  and  triple  scattering  for  2-D  and 
3-0  media  in  the  case  of  isotropic  scattering.  Gao  et  al.  (1983,  1984) 
derived  up  to  seventh  order  scattering  and  then  obtained  the  approximate 
formulas  of  multiple  scattering  in  time  domain  for  2-D  and  3-D  media  using 
curve-fitting  technique.  However,  the  formulas  derived  are  for  the  case  in 
which  the  source  and  sensor  are  located  in  the  same  point.  On  the  other  hand, 
the  most  proiainent  evidences  of  aiultiple  scattering  would  be  manifested  if  the 
sensor  could  be  situated  at  some  place  between  the  source  and  the  point  apart 
from  the  source  by  one  mean- free -path  of  scattering  (this  will  be  shown 
later).  Therefore  it  may  be  difficult  to  use  these  formulas  for 
discriminating  the  scattering  attenuation  from  the  intrinsic  attenuation, 
though  the  formulation  may  be  very  useful  in  other  calculations. 

In  this  paper,  we  derive  the  formulation  of  seismic  energy 
transfer  under  multiple  scattering  by  using  the  radiative  transfer  equation 
technique  developed  in  the  aatrophysical  optics  and  the  neutron  transport 
theory  and  explore  the  possibilities  of  using  this  approach  to  separate  the 
scattering  and  intrinsic  attenuation. 

Historically,  multiple  scattering  theory  has  been  developed  along  two 
independent  approaches:  the  analytic  theory  and  the  transport  theory  (for 
review  see  Ishlmaru  1977).  Both  are  based  on  the  statistical  treatment  of 
wave  propagation  in  random  media.  Because  the  complex  heterogeneities  are 
■odeled  with  a  random  medium,  the  waveflelds  propagating  therein  are  also 
random  waveflelds.  He  are  interested  only  in  some  statistical  quantities  of 
the  wavefield,  such  as  the  mean  intensity,  phase  and  amplitude  fluctuations, 


various  correlation  functions,  pulsa  spreading,  angular  broadening,  etc.  All 
of  these  quantities  can  be  obtained  from  the  moments  of  the  random  field.  The 
analytic  theory  starts  with  basic  differential  equations  such  as  wave 
equations  and,  by  introducing  the  scattering  and  absorption  characteristics  of 
the  random  heterogeneities,  derives  the  differential  or  integro-differential 
equations  for  the  moments  of  the  wavefields.  There  are  basically  two  branches 
in  the  analytic  theory:  the  renormalization  method  and  the  small-angle 
approximation  method.  In  the  first  branch  the  renormalization  procedure  was 
used  for  tne  formal  perturbation  series  and  the  exact  equation  for  the  first 
moment  (the  mean  field),  known  as  the  Dyson  equation,  and  for  the  second 
moment  (the  correlation  function),  the  Bethe-Salpeter  equations  were  derived. 
These  equations  are  exact  in  the  sense  that  the  multiple  scattering  of  all 
orders,  as  well  as  the  diffraction  and  interference  effect  are  all  included  in 
the  equations.  However,  since  the  operator  involved  in  these  equations  are  in 
the  form  of  infinite  series,  there  is  no  solution  aveilable  at  present. 
Approxiaa tions  have  to  be  made  to  the  operator  before  some  practical  solutions 
can  be  obtained.  The  most  widely  used  approximation  is  the  first  order 
smoothing  approximation  as  called  by  Frisch  (1968)  (see  also  Ishimaru  1978,  v. 
2),  in  which  the  local  Born  approximation  of  the  fluctuating  field  (or 
equivalently  the  Bilocal  approximation  to  the  mean  field)  is  applied  to  the 
Dyson  equation  and  the  ladder  approximation  is  applied  to  the  Bethe-Salpeter 
equation.  These  approximations  can  be  obtained  by  either  the  Feynman  diagram 
method  or  the  Bogoliubov  smoothing  method  in  the  operator  form  (Frish  1968, 
Tatarskii  1971,  Ishimaru  1978;  for  the  various  names  of  the  first  order 
smoothing  approximation,  see  also  Wu  1982b,  footnote  2).  The  justification 
for  the  use  of  this  approximation  has  been  clarified  by  Frisch  (1968)  by 
introducing  the  generalized  Reynolds  number.  The  basic  physical  condition  for 


the  valid  use  of  the  approximation  is  the  scattered  field  within  a  correlation 
length  being  weak  compared  with  the  incident  field.  In  the  case  of  large 
scale  inhomogeneities i  Fante  (1982)  has  shown  that  a  sufficient  condition  for 
applying  the  ladder  approximation  is  the  mean  free  path  for  multiple 
scattering  being  large  in  comparison  with  the  correlation  length  of  the 
medium.  This  condition  is  usually  satisfied  in  the  context  of  seismic  wave 
scattering  in  the  lithosphere.  The  first  order  smoothing  approximation  to  the 
Dyson  equation  and  Bethe-Salpeter  equation  can  be  shown  (Frisch  1968)  to  be 
equivalent  to  the  Fold*— Twersky  system  of  equations,  which  have  been  developed 
independently  for  discrete  random  media,  i.e.  the  media  with  randomly 
distributed  scatterers.  There  are  still  no  general  solutions  for  these 
equations  and  further  approximations  are  needed  to  put  them  into  practi  .1 
use.  For  small  size  inhomogeneities,  there  are  some  general  solutions  for  cne 
mean  field,  but  nb  useful  results  for  the  second  moments  (Tatarski  1971,  §61, 
Ishimaru  1978,  ch.  14).  It  has  been  shown  that  the  first  order  smoothing 
approximation  of  the  Dyson  and  Bethe-Sapeter  equations  can  lead  to  a  radiative 
transfer  equation  for  the  specific  intensity  which  is  the  3D  spatial  Fourier 
transform  of  the  spatial  correlation  function  of  the  wavefield  when  the 
correlation  function  is  a  slowly  varying  function  in  space  (Barabanenkov  1969, 
1971,  Tatarskii  1971,  §63,  Ishimaru  1975,  1978).  Similarly,  a  generalized 
radiative  transfer  equation  can  be  derived  for  the  frequency  correlation 
function  (Ishimaru  1978).  Thereby  the  link  bridge  has  been  established  between 
the  analytic  theory  and  the  transport  theory. 

The  second  branch  of  the  analytic  theory  includes  all  the  small-angle- 
scattering  methods.  Because  of  the  small  scattering  angle  approximation  or 
forward-scattering  approximation,  the  basic  starting  point  of  the  method  is 
the  parabolic  wave  equation.  There  are  two  approaches:  parabolic  equation 


approach  and  Feynman  path  integral  approach.  Tatarskii  applies  the  Markov 
approximation  to  the  parabolic  wave  equation,  so  the  theory  of  Markov  process 
can  be  used  to  the  study  of  the  problem  (Tatarskii  1971).  Uscinski,  on  the 
other  hand,  uses  the  plane  wave  decomposition  and  phase-screen  technique  to 
the  parabolic  wave  equation  (Uscinski  1977).  At  present,  the  parabolic 
equation  methods  can  have  only  approximate  solutions  for  up  to  the  fourth 
moment  equations.  The  path-integral  approach  starts  with  the  Feynman 
path-integral  representation  of  the  parabolic  wave  equation  and  makes  use  of 
the  small  scattering-angle  approximation  and  Markov  approximation  (Dashen 
1977,  Platte  et  al. ,  1979).  It  can  obtain  solutions  for  any  higher  order 
moments  for  the  Gaussian  statistics.  Flatte  et  al.  have  applied  this  approach 
to  the  ocean  acoustics  and  obtained  the  expressions  for  phase  and  intensity 
fluctuations,  various  correlations  and  pulse  wandering  and  spreading  etc. 

The  transport  theory  (or  radiative  transfer  theory)  is  a  phenomenological 
approach.  It  does  not  start  with  the  wave  equation,  but  deals  directly  with 
the  energy  transport  process.  Therefore,  only  energy  or  intensity  arithmetic 
appears  in  the  theory  and  no  wave  interference  is  considered.  This  treatment 
much  simplifies  the  mathematics.  Historically  it  appeared  earlier  than  the 
analytic  theory,  and  has  its  root  from  Boltzmann's  equations  in  the  kinetic 
theory  of  gases  and  in  the  neutron  transport  theory.  It  was  introduced  into 
astrophysical  optics  by  Schuster  (1905),  Chandrasekhar  (1950)  and  others  and 
is  now  widely  used  in  the  multiple  scattering  treatment  in  the  astrophysical 
optics,  ocean  acoustics,  neutron  transport  theory,  electromagnetic  wave  remote 
sensing,  marine  biology,  etc.  (Chandrasekhar  1950,  Sobolev  1963,  Menzel  1966, 
Davison  1958,  Bell  and  Glasstone  1970,  Flatte  1979,  Kong  et  al.  1984, 

Jerlov  1976).  This  approach  also  has  its  shortcomings.  It  can  only  deal  with 
the  second  moments,  it  does  not  account  for  the  diffraction  and  interference 


phenomena.  However,  there  are  some  new  developments  recently,  which 
incorporate  some  wave  interference  effects  into  the  radiative  transfer 
equation.  For  example,  in  deriving  the  transfer  equations  from  the 
Bethe-Salpeter  equation,  beside  the  ladder  terms  (which  alone  will  lead  to  the 
regular  intensity  transfer  equation),  the  cyclical  diagrams  are  also  included, 
resulting  in  a  modified  radiative  transfer  equation,  which  can  account  for  the 
backscattering  enhancement  due  to  the  constructive  interference  effect  caused 
by  the  double  passage  of  the  backscattered  waves  (Zuniga  et  al.  1980). 
So-called  "wave  radiative  transfer  theory"  based  on  the  second  order 
approximations  to  the  Bethe-Salpeter  equation  is  also  under  development  (Tsang 
and  Ishimaru  1983). 

For  the  coda  envelopes  or  coda  energy  problems  of  local  earthquakes,  it 
is  apparently  a  wide-angle  scattering  problem,  so  that  the  transport  theory  is 
probably  the  most  effective  method  to  treat  it  at  present.  In  this  paper  we 
use  the  frequency  domain  formulation  mainly  from  the  neutron  transport  theory 
and  the  electromagnetic  wave  propagation  (Davison  1958,  Liu  and  Ishimaru  1974, 
Fante  1973,  Ishimaru  1978)  to  the  energy  density  decay  with  distance  of  the 
seismic  waves  from  local  earthquakes,  and  discuss  the  possibility  of  using  the 
decay  curves  to  evaluate  the  relative  strengths  of  the  intrinsic  absorption 
and  the  scattering  coefficient  of  the  medium  in  the  region  studied.  Some 
examples  are  given  for  the  Hindu-Kush  region.  The  results  and  their 
geophysical  meaning  are  also  discussed. 

2.  DEFINITIONS  AND  NOTATIONS 

It  is  difficult  to  keep  all  the  notations  and  terminology  in  radiative 
transfer  theory  without  causing  ambiguities  and  contradictions  with  the 
traditional  notation  and  terminology  in  seismology,  when  the  theory  is 
introduced  into  seismology.  I  will  basically  follow  Ishimaru  (1978)  and  make 
some  necessary  changes  to  keep  the  notations  self-consistent. 


I(jr,Q):  Specific  intensity  or  directional  Intensity.  It  Is  the  most 


fundamental  quantity  In  transport  theory.  It  gives  the  power  flowing 

A  A 

within  a  unit  solid  angle  in  the  direction  Q,  here  Q  is  the  unit  vector, 

A 

emanated  from  a  unti  area  perpendicular  to  Q,  in  a  unit  frequency  band. 
The  specific  intensity  is  defined  for  a  fraquency  o>,  which  is  omitted  in 
the  notation. 

In  this  paper  we  consider  the  S  wave  and  its  coda  for  small  local 
earthquakes.  Since  the  P  wave  energy  is  much  smaller  than  the  S  wave 
energy  for  a  double-couple  point  source  which  is  the  source  model  for 

A 

small  earthquakes,  we  consider  here  l(jr,Q)  as  only  the  S  wave  energy  by 
neglecting  the  mode  converted  energy  from  P  waves.  We  assume  here  also 

A 

that  the  wave  energy  described  by  I(jr,Q)  is  depolarized,  i.e.  the 
energy  is  equally  partitioned  between  the  two  orthogonal  components  of  S 
waves.  This  agrees  generally  with  the  observations.  Because  of  the  free 
surface  reflection  and  the  scattering  by  heterogeneities,  the  S  waves 
from  a  double-couple  source  get  quickly  depolarized.  From  the  results  of 
this  paper,  the  energy  density  decay  curves  for  the  two  orthogonal 
components  are  very  similar  to  each  other,  which  further  validate  the 
aasump  tlons . 

In  order  to  measure  the  specific  intensity  (or  directional  intensity), 
we  need  strongly  directional  sensors,  which  are  not  available  in  the 
seismological  practice.  Therefore  the  specific  intensity  is  not  the 
quantity  measured  in  practice,  but  is  the  Important  concept  and  quantity 
for  theoretical  derivations. 

I(£)s  Average  intensity,  defined  by 

"  h  L  I(r'°)dQ’ 

is  the  intensity  at  point  £  averaged  over  all  directions. 


E(r):  Energy  density,  defined  by 

E(r)-i/  I(r,Q)dQ  -  i<r)  .  (2.2) 

—  4ic  —  — 

where  C  is  the  wave  velocity. 

^(r):  Flux  density  vector,  defined  by 

A  A 

J(r)  -  /  I(r,Q)QdQ  .  (2.3) 

4k 

A  A 

The  net  flux  density  in  a  particular  direction  Q0  is  defined  as  Q0.J(r). 

A 

It  is  the  net  power  transferred  along  the  Q0  direction  across  a  unit 

A 

area  perpendicular  to  Q0.  In  this  paper,  we  also  use  the  notation  for  the 
energy  flux  density,  i.e.  the  power  flux  density  divided  by  the  wave 
velocity  c. 

A  A 

S(Q,Q0):  Scattering  intensity  function  of  a  random  medium,  which  is  related 

A  A 

to  the  single  scattering  amplitude  f(0,Qo)  of  an  elementary  volume  dV  of 
the  inhomogeneou.s  medium  by 

<|f(Q,Q0)|2> 

S(fl,Q0)  - - — - ,  (2.4) 

A  A 

where  <  >  denotes  taking  ensemble  average.  S(Q,Q0)  gives  the  scattered 

A 

power  in  Q  direction  within  a  unit  solid  angle  by  a  unit  volume  of  the 

A 

random  medium  for  a  unit  flux  density  of  incident  wave  in  Q0  direction. 

In  this  paper  we  will  give  a  unified  treatment  for  both  the  discrete 
and  the  continuous  random  media.  For  a  discrete  random  medium  composed 

A  A 

of  randomly  distributed  scatterers,  S(Q,Q0)  is  defined  by  the  scattering 
characteristics  of  individual  scatterers;  while  in  the  case  of  random 
continua,  we  can  choose  the  volume  elements  small  enough  so  that  we  can 

A  A 

derive  the  single  scattering  amplitude  f(Q,Q0)  by  the  Born  approximation. 


g(Q,Q0):  Directional  scattering  coefficient,  defined  by 

A  A  A  A 

g<0,flo)  -  4itS(Q,Q0)  .  (2.5) 

for  the  definition  and  the  derivation  for  elastic  random  media,  see 
paper  II  (Wu  and  Aki,  1984b). 
r|a=g:  Scattering  coefficient  of  the  medium  defined  by 

Hs  “  /  S(Q,Q0)dQ  -j-f  g(0,Qo)dQ  ,  (2.6) 

4*  4k 

which  gives  the  total  power  loss  due  to  scattering  by  a  unit  volume 
random  medium  per  unit  flux  density  of  Incident  wave  under  the  single 
scattering  assumption. 

TJs=b:  Absorption  coefficient  of  the  medium,  which  gives  the  power  loss  due 
absorption  by  a  unit  volume  random  medium  per  unit  flux  density  of 
incident  wave. 

T)as  Extinction  coefficient  of  the  medium,  defined  by 

T}*-na+ns  (2.7) 

-i0=a:  Correlation  length  of  the  random  medium. 

I*a  ■  1/t|«:  Extinction  length  of  the  medium. 

I*a  "  1/Da8  Absorption  length  of  the  medium.  (2.8) 

L®  ■  1/t)s8  Scattering  length  or  scattering  mean  free  path  of  the  medium. 

D®s  Numerical  extinction  distance,  which  is  called  "optical  distance"  in  optics. 
D®s  Numerical  absorption  distance, 

D®:  Numerical  scattering  distance,  defined  by 
D®  -  r/Le, 

D®  ■  r/La,  (2.9) 

D,  -  r/Lg, 


where  r  is  the  travel  distance 


B0:  Medium  seismic  albedo,  defined  by 


*ls 


B°  "  "  na+*la 


(2.10 


D(Q,Q0):  Scattering  directivity,  defined  by 

A  A  A  A 

g(Q,Q0)  4itS(0,Qo) 


A  A 

D(Q,Q0)  - 


(2.11. 


U*  Tie 

It  Is  the  normalized  directional  scattering  coefficient,  and  satisfies 


/  d<q.q0)  -  i. 


(2.12) 


4it 


that  means  its  average  over  all  the  directions  is  equal  to  unit.  In 
the  case  of  isotropic  scattering 

A  A 

D(Q,Q0  5  1  .  (2.13) 

Its  relation  with  the  "phase  function"  in  the  radiative  transfer  theory 
(Chandrasekhar  1930',  Ishinara  1978)  is 

A  A  A  A 

D(Q,Q0)  -  BoP (fi,Q0)  .  (2.14, 

A  A 

p(Q,Q0):  Phase  function  (see  2.14). 

In  the  case  of  a  discrete  random  medium  having  statistically  uniformly 
distributed  random  scatterers  with  number  density  n,  we  have 

A  A 

0d(Q,Qo):  Differential  (or  directional)  scattering  cross-section  of  the 
scatterers. 

A  A  A  A 

S(Q,Q0)  -  nod(Q,Q0).  (2.15) 

as:  Scattering  cross-section  of  the  scatterers.  defined  by 

(2.16) 


Os  m  I  0d(Q,Qo)dQ 

4n 

oa s  Absorption  cross-section  of  the  scatterers. 


To tal  cross-section  of  the  scatterers 


3.  ENERGY  DENSITY  DISTRIBUTION  IN  THE  CASE  OF  ISOTROPIC  SCATTERING 


Knowing  the  extinction  coefficient  and  scattering  coefficient  of  the 

A  A 

aedlua  t)s  and  the  scattering  directivity  D(Q,Q0)  or  the  scattering 

A  A 

intensity  function  of  the  medium  S(Q,Q0)  defined  by  (2.6),  (2.7),  (2.11)  and 
(2.4),  we  can  obtain  the  differential  equation  for  the  specific  intensity 

A 

JC(r,Q),  the  "equation  of  transfer"  (Chandrasekhar  1950,  I,  Ishimara  1978,  ch. 

7)  : 

A 

dI(r,Q)  A  A  A 

— 3? -  -  -  Tie*<£>a>  +  /  S(Q,0o)I(r,Qo)dQo  +  W(r,Q) 

4 x  —  — 

-  -  t1aI(r,Q)  +  ^  /  D(Q,Q0)I(r,Q0)dQ0+W(r,Q)  ,  (3.1) 

—  4* 

A 

where  V(x_,Q)  is  the  source  intensity  function,  which  defines  the  amount  of 

A 

power  emitted  from  the  sources  into  the  direction  Q  per  unit  solid  angle.  In 
(3.1),  di  is  the  length  of  a  cylindrical  elementary  volume  of  unit  cross 

A 

section  in  the  medium  with  the  axis  of  the  cylinder  in  Q  direction  (Fig.  3.1). 
Therefore  the  left  hand  side  of  (3.1)  represents  the  total  change  of  the 
specific  intensity  for  a  unit  travel  distance.  The  first  term  in  right  hand 

A 

side  of  (3.1)  is  the  loss  of  power  in  Q  direction  due  to  absorption  and 
scattering,  whereas  the  second  term  gives  the  gain  of  power  in  that  direction 
from  the  scattered  waves  for  the  Incident  Intensity  from  all  directions  and 
the  third  term  is  the  energy  supply  from  the  sources.  No  general  analytic 
solutions  are  available  for  (3.1).  Some  methods  such  as  the  Gauss-quadrature 
can  be  used  to  obtain  the  numerical  solutions  for  a  general  scattering 
function.  Let  us  first  consider  the  simplest  case  of  isotropic  scattering. 

A  A 

In  this  case  the  scattering  directivity  D(Q,Q0)=1.  Integrating  (3.1)  over  all 

A 

directions  Q,  we  obtain  equation  for  the  average  intensity  l(r)  or  the 


energy  density  E(r)  (2.2) 

dE(r)  ,  t>g 

“3J —  -  -  TieE(r)+£  /  /  I(r,Q0)dQ0  +  W(r,Q)]dQ 

—  4s  4s  —  — 

-  -  q€E(r)  +  Q(A),  (3.2) 

where  C  is  the  wave  velocity.  (3.2)  is  in  a  form  of  first  order  differential 

equation,  in  which  the  second  tern  in  RHS  is  the  source  term 

QU)-4j  [— J  I(r,Q0)dQ0  +  w(r,Q)]dQ  .  (3.3) 

4s  4s  —  — 

The  general  solution  for  (3.2)  is 

X 

E(r)  -  Ae"V  +  /  QUi)e-ne(A-A1)d2L  ,  (3.4) 

o 

where  A  is  a  constant. 

The  energy  density  (3.4)  is  composed  of  two  terms.  The  first  term  is  a 
simple  exponential  decay  with  the  extinction  coefficient  T)e  as  its  attenuation 
coefficient;  this  is  the  coherent  energy  density  Ec  or  "reduced  energy 
density"  (Ishimaru  1978).  The  second  term  is  therefore  the  diffuse  energy 
density  Eg  which  is  produced  by  scattering.  Applying  the  initial  condition 

E(£0)  -  Eln,  (3.5) 

where  E^n  is  the  Incident  energy,  we  get 

E(r)  -  Ec(r)  +  Ed(r) 

EcOr)  •  Ein  e^e* 

1 

Eg(r)  -  /  Q(Jl|,)e_T>e(*"*i)  dJlj 

/  /  ItJ  /  I(r ,Q0)dQ0  +  W(r,Q)]e“Tle<*-VdQd*l  . 

o  4s  4s  -  - 


l 

C 


(3.8) 


Ia  order  to  calculate  the  diffuse  term  (3.8),  we  need  to  know  the  intensity 

A 

l(r»Q0)  which  is  related  to  the  total  energy  density.  Therefore  (3.8)  is  in 
the  form  of  integral  equation.  To  carry  out  the  integration  with  respect  to 

A  A 

Q,  we  note  that,  the  intensity  gain  in  the  direction  Q  within  dQ  are 
contributed  from  the  intensity  of  all  the  volume  elements  dV^  at  _£]_  within  the 
elementary  solid  angle,  and 

dVi  -  dQlr-rj j2dli  .  (3.9) 

Therefore  (3.8)  becomes 

,  .  a.  A  e"T'elr“rll 

Ed(r)  -  /  [t|aE(r1)+—  W(rltQ)]  - —  dVL  .  (3.10) 

-  V  -  c  -  Awlr-ril2 

The  integration  is  over  the  volume  of  the  random  medium.  The  integral 
equation  for  the  total  energy  density  becomes  (see  also  Ishimaru,  1978,  ch. 
12). 

•  A 

E(r)  -  2i„e-M  +  /  Ins2(rl)4€(ri,fi))G0(r-r1)dV1  ,  (3.11) 


where 


4x 


e(r,Q)  -~V  (ri,Q) 

is  the  source  energy  density  function,  and 


(3.12) 


G0(r-ri)  -T 


■T|,R  r*ri ' 


4itR<^  .  |  ,2 

4w|r-ri| 


(3.13) 


Integral  equation  (3.11)  can  also  be  derived  from  the  first  order 
smoothing  approximation  of  the  Dyson  and  Bethe-Salpeter  equations  (Lin  and 
Ishimaru  1974). 

From  (3.11),  the  energy  density  E(r)  is  totally  defined  by  the  Incident 


field,  the  source-function,  and  the  volume  of  the  random  medium.  For  the 
problems  of  seismic  coda  waves  of  local  earthquakes,  the  distances  between  the 
stations  and  the  sources  are  short  compared  to  the  travel  time  of  coda  waves. 


21. 


As  the  first  approximation,  we  consider  the  problem  of  a  point  source  Located 
in  an  infinite  randomly  inhomogeneous  medium.  The  effect  of  the  free  surface 
is  like  a  mirror  reflecting  the  half  random  space  to  a  whole  random  space  with 
the  upper  half  space  being  the  mirror  image  of  the  lower  half  space.  The 
limited  thickness  of  the  Lithosphere,  which  is  supposed  to  be  more 
heterogeneous  than  the  asthenosphere  beneath  will  have  influence  on  the  coda 
of  Later  part.  Further  discussion  about  the  Limitation  of  the  model  will  be 
given  Later  in  this  paper. 

In  (3.LL),  suppose  the  incident  field  E^n  ■  0  and  the  point  source  is 
Located  at  _r  ■  0,  radiating  the  total  power  P0.  Then 

*o 

e(r)  -  —  6(r)  -  E06(r)  (3.14) 

The  equation  (3.LL)  becomes 


E<rJ  *  Eo  +  /  E(rl> 


e-He | r-ri j 


“  4*  tern 


dVl 


-  E0  G0(r)  +  /  risE(£1)G0(£-£i)dV1  .  (3.15) 

V 

This  is  a  Faitung  type  or  -convolution  type  integral  equation 

(Trieomi  1957,  Carrier  et  al.  1966),  Fourier  transform  method  can  be  used  for 

solution.  Assuming  E0  ■  1,  the  solution  can  be  written  as  (see 

Davison  1958,  Lin  and  Ishimaru  1974,  Ishimaru  1978  (12-21)) 

T)€Pd  n«  - 

E'r}  “  4 TT  exP<-T1edor)  +  4^7  /  f(s,B0)exp(-ners)ds 


Ed(r)  +  Ec(r) 


(3.16) 


where 


2  d20(l-d20) 

Pd  - - - -  ,  (3.17) 

Bo<<»  o+Bq-D 

end  dQ  is  the  diffuse  multiplier  determined  by 

B0  l+d0 

2d7An  (i^;) "  1{  (3-l8> 

end 

B  B 

f(s,B0)  -  {[1-  ~  tenh-l(i)]2  +  (7 -AV1  .  (3.19) 

S  S  AS 

The  first  term  in  (3.16)  is  the  diffuse  term  £4,  which  is  attributed  to  the 
pole  residue  in  the  complex  spatial  frequency  plane,  and  the  second  term, 
coherent  term  Ee  is  from  the  branch  cut  integration. 

Fig.  3.2  shows  the  relation  between  the  diffuse  multiplier  dQ  and  the 
medium  albedo  B0.  dQ  is  always  less  than  1.  When  distance  r  is  large, 
especially  for  large  Bot  the  diffuse  term  becomes  dominant  (see  also  Fig.  9) , 
and  E(r)  will  be  approximately  an  exponential  decay  with  an  apparent 
attenuation  coefficient  d0T]e,  which  is  less  than  the  extinction  coefficient 
t|€.  The  degree  of  reduction  depends  on  the  albedo  BQ.  The  diffuse  terra  can 
also  be  written  as 

’le^d 

Ed<r>  "  4^7-  exp[-(ria+dsris)rj, 
d0-(i-B0) 

<*s  -  - S -  .  (3.20) 

ao 

ds  is  a  multiplier  and  d8n8  gives  the  effective  contributions  of  the 
scattering  coefficient  to  the  apparent  attenuations.  ds  is  also  plotted  in 
Fig.  2.  Table  1  lists  some  values  of  d0  and  ds  versus  B0. 


23. 


The  coherent  tern  can  also  be  written  as 

Tie  1  T)er  at 

Ec<r)  -4^/  f^.BoJexpC-r-)  Sk  ,  (3.21) 

o  *  5 

by  setting  £  “  1/s  for  the  convenience  of  computation.  Fig.  3.3  shows  the 
behavior  of  the  two  factors  of  the  integrand  for  different  numerical 
extinction  distances  Da  -  ri€r  and  different  medium  albedo  B0.  exp(-De/£)/£2 

has  a  sharp  peak  for  small  De  when  £  is  small;  whereas  f(£,B0)  is  nearly 
singular  for  small  BQ  when  £  1-s  close  to  1.  Therefore,  in  doing  numerical 
integration,  we  used  Romberger  integration  method  for  three  separate  segments 
to  take  care  of  the  abrupt  changes  of  the  integrand  at  both  ends  of  the 
interval.  The  Gauss-Legendre  quadrature  is  also  used  to  check  the  results. 

It  turned  out  that  the  Gauss-Legendre  quadrature  of  order  10  gives  fairly  good 
results. 

In  the  following  we  will  show  some  numerical  results  of  the  energy 
density  distribution  along  the  travel  path  from  the  source  point.  In  the  case 
of  homogeneous  media,  the  decay  of  energy  density  with  distance  is  only  due  to 
geometric  spreading.  For  a  isotropic  point  source,  the  decay  is  l/4xr2. 
Therefore,  we  normalize  the  distribution  for  inhomogeneous  media  (3.16)  by  the 
homogeneous  distribution,  i.e.  multiply  both  sides  of  (3.16)  with  4xr2, 

I  jp 

En(r)  -  4wr2E(r)  «  ilePdr  exp(-ried0r)+ner/f (£,B0)exp(— r-)^7 

o  ^  4 

-  DePdexp(-d0De)+De  /  f (£,B0)exp(-De/£>^-  , 

o  £2 


(3.22) 


where  E^r)  stands  for  the  normalized  energy  density  distribution.  Fig.  3.4 
gives  the  results  for  different  medium  albedo  B0.  The  diffuse  term  and  the 
coherent  term  are  also  plotted  in  the  figure  for  comparison.  The  coherent 
term  has  little  changes  for  different  B0,  whereas  the  diffuse  term  varies 
dramatically  with  BQ,  especially  when  Bo>0.5,  i.e.  when  scattering  is 
dominant.  This  gives  the  possibility  of  using  the  energy  density  decay  curves 
to  calculate  the  extinction  coefficient  rie  and  the  medium  albedo  B0,  hence  to 
separate  the  absorption  coefficient  T)a  and  the  scattering  coefficient  T)s.  In 
the  case  of  Bo>0.5,  the  diffuse  term  is  dominant.  There  will  be  a  peak  on  the 
E(r)  curve,  the  position  of  the  peak  will  depend  on  t)e  and  B0  of  the  medium. 
When  Bo<0.5,  the  coherent  term  is  dominant  for  De<2.  Therefore  the  shape  of 
the  curve  is  not  very  sensitive  to  the  change  of  B0,  so  that  the  separation  of 
scattering  from  absorption  becomes  difficult. 

By  assuming  a  point  source  with  Eg"!,  we  get  E(r)  around  the  peak  with 
values  greater  than  1,  that  need  some  explanation.  As  shown  in  Fig.  3.5,  the 
normalized  energy  density  E,j(r)  ■  4«r2E(r)  represents  the  energy  received  by 
the  ring  shell  (hatched).  In  a  homogeneous  medium,  if  there  is  no  absorption, 
the  energy  received  will  be  equal  to  the  source  energy.  In  a  scattering 
medium,  the  wave  energy  can  go  outward  and  inward  across  the  shell.  We  denote 
the  outward  energy  flux  by  Fr+  and  the  inward  energy  flux  by  Fr~.  In  the 
figure,  we  sketched  one  possible  path  of  multiple  scattering.  No  matter  how 
complicated  the  path  is  and  how  long  the  time  delay  is  compared  to  the  direct 
path,  the  closed  ring  shell  will  eventually  receive  all  the  energy  emitted  by 
the  source.  There  is  no  escape!  Therefore,  in  this  case  the  Fr+  is  equal  to 
the  total  energy.  However,  the  shell  will  also  receive  the  inward  scattered 
energy,  so  the  total  received  energy  Fr++Fr"is  greater  than  E0.  Of  course  the 
net  energy  flux  Fr+-Fr“is  always  less  than  E0. 


If  there  exists  absorption 


Che  amount  of  received  energy  will  depend  on  the  energy  balance  between  the 
absorption  loss  and  the  inward-scattering  gain.  Near  the  source,  r  is  small, 
the  ring  shell  has  a  small  surface  area  for  receiving  the  inward-scattered 
energy,  so  E^r)  *  E0.  When  r  increases,  the  surface  area  of  the  shell  also 
increases,  so  that  more  inward-scattered  energy  can  be  received,  resulting  in 
the  growth  of  E^r).  However,  the  absorption  loss  also  grows  with  r  due  to 
the  increase  of  the  path  length.  Upon  to  some  distance  r,  the  growth  rate  of 
gain  is  equal  to  the  growth  rate  of  loss,  the  curve  reaches  its  maximum. 

Beyond  this  distance,  the  absorption  loss  prevails. 

Fig.  3.6  replots  the  curves  of  Fig.  4  in  a  semi-logarithm  coordinate 
system.  Fig.  7  and  8  plot  some  En(r)  curves  for  cases  of  constant  absorption 
and  constant  scattering  respectively.  In  this  paper  b=T)g  g=T)s.  Fig.  3.7 
shows  the  influence  of  different  scattering  coefficients  on  the  energy  density 
distribution  curve*  of  a  constant  absorption  medium.  The  distance  is 
normalized  by  the  absorption  length  of  the  medium  La  ■  l/-pa.  It  is  seen  from 
the  figure  that,  for  large  distances  compared  with  the  absorption  length  of 
the  medium,  the  decay  of  the  energy  density  is  nearly  exponential  with  an 
apparent  attenuation  coefficient  different  from  both  the  extinctions 
coefficient  and  the  absorption  coefficient.  In  the  figure,  b  is  the  true 

A 

absorption  coefficient,  b  is  the  apparent  attenuation  coefficient  measured 
from  the  slope  of  the  curve.  It  can  be  seen  that,  for  strong  scattering 
(Bo>0.5),  the  apparent  attenuation  is  much  bigger  than  the  absorption 

A 

coefficient  but  much  smaller  than  the  extinction  coefficient  (for  B0  =  0.9,  b 
■  4.5b  ”  0.45  T)e).  For  weak  scattering  (Bo<0.5),  the  influence  of  scattering 
to  the  apparent  attenuation  is  less  appreciable.  When  B0  -  0.5,  b  =  1.62b. 

On  the  other  hand,  for  small  absorption  distance  (Da<l),  the  shape  of  E(r) 
curve  varies  drastically  depending  on  the  values  of  B0,  which  provides  the 
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basis  for  the  separation  of  scattering  effect  and  the  absorption  effect.  Fig. 
3.8  in  a  similar  way  shows  the  influence  of  absorption  on  the  E(r)  curve  of  a 
constant  scattering  medium. 

In  order  to  compare  the  relative  contributions  of  the  diffuse  term  and 
the  coherent  term,  we  plot  them  on  Fig.  3.9  and  Fig.  3.10  with  the  distance 
normalized  by  the  extinction  length  Le  and  scattering  length  Ls  respectively. 

Now,  we  will  derive  the  radial  energy  flux  density  Jr(r).  We  know  the 


energy  conservation  relation  (see  Ishimaru  1978,  (7.28)) 


"a  ,  »  i  ,  * 

div  J(r)  -  -  7-  /  I(r,Q)dQ  +ff  W(r,Q)dQ, 


(3.23) 


4it  — 


Ait  — 


where  _J(r)  is  the  energy  flux  density  vector,  C  is  the  wave  velocity  and 

A 

W(r,Q)  is  the  source  intensity.  For  isotopic  scattering  in  the  source  free 


region 


div  J[(jr)  »T}aE(r).  (3.24) 

In  view  of  the  spherical  symmetry,  there  is  no  transverse  component  of  J(r) , 


therefore  (3.24)  becomes 


J(r)  -  (r2Jr)  -  -  riaE(r)  . 


J r  “  "  7T  /  E(r)r2dr  -  -p-  j  E(r)r2dr 


(3.25) 


(3.26) 


Normalizing  Jr  by  the  homogeneous  case,  we  get 


Jnr(r)  “  4xr2Jr(r)  -  T}a  /  4nr2E(r)dr  «  t ia  J  En(r)dr  .  (3.27) 


«*-  ■**-  ■*  -  .*  *  .*  *-*  '  -  ■  '  ■  "J 


.  '  .1 
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Substituting  (3.22)  into  (3.27)  yields 
Jnr(r)  *  4wr2Jr(r) 

Pd  ,  1  °e 

-  (l-B0){T-(De+i-)exp(d0De)  +  /  f(5,B0)(l-  r-)exp(-De/5)dS  (3.28) 

a°  O  6 

Fig.  3.11  and  3.12  give  some  numerical  results  with  the  distance  normalized  by 
the  extinction  length  and  by  the  absorption  length  respectively,  together  with 
the  results  for  the  forward  scattering  approximation  (see  next  section).  It 
can  be  seen  that  the  radial  net  flux  is  always  smaller  than  the  source  energy 
E0.  However,  the  radial  energy  flux  is  difficult  to  measure  in  the  practice 
of  seismology.  The  reason  is  the  difficulty  of  separating  the  inward  and 
outward  energy  flow.  Nevertheless,  the  comparison  between  E(r)  and  Jr(r) 
helps  us  understand  the  multiple  scattering  process. 


4.  STRONG  FORWARD  SCATTERING:  THE  CASE  OF  LARGE  SCALE  INHOMOGENEITIES 

From  the  analysis  of  coda  generations  for  local  earthquakes,  we  conclude 
that  the  lithosphere  in  tectonically  active  regions  may  be  rich  in  small  scale 
heterogeneities  (less  than  1  km)  (paper  II).  On  the  other  hand,  by  measuring 
the  phase  and  amplitude  fluctuations  in  large  seismic  arrays  as  LASA  and 
NORSAR,  large  scale  velocity  inhomogeneities  (10-20  km)  underneath  the  arrays 
were  revealed  (Aki  1973,  Capon  1974,  Berteusson  et  al.  1975).  Therefore,  the 
lithosphere  may  have  multi-scale  inhomogeneities.  For  short  period  seismic 
waves  (around  1  Hz),  the  scattering  by  the  small  scale  heterogeneities  may  be 
in  the  Rayleigh  and  Mie  scattering  region.  From  the  elastic  scattering 
pattern  (paper  I,  II),  we  may  approximately  use  the  isotropic  scattering 
approximation.  However,  for  the  large  scale  velocity  inhomogeneities,  the 
forward  scattering  is  dominant.  The  energy  density  distribution  with  distance 
will  be  quite  different  from  the  case  of  isotropic  scattering.  Since  most  of 
the  scattered  energy  is  concentrated  in  the  forward  direction  within  a  small 
cone,  the  focussing  and  defocussing,  diffraction  interference  phenomena  become 
important.  Most  of  the  scattered  energy  arrives  at  the  receiver  point  with 
much  shorter  travel  paths,  so  that  the  energy  delay  due  to  scattering  is  much 
less  severe  than  the  case  of  isotropic  scattering.  From  a  reasoning  similar 
to  that  in  Fig.  3.5,  we  can  see  that,  the  normalized  energy  density  decay 
curve  will  not  have  a  peak  of  value  greater  than  1.  Because  the  inward 
scattered  energy  is  much  less  than  the  outward  scattered  energy,  the  energy 
density  which  is  Jr+  +  Jx*,  where  Jr+  and  Jr_  is  the  outward  and  inward  radial 
energy  flux  respectively,  will  have  no  too  much  difference  from  the  net 
energy  flux  Jr  ■  Jr+-Jr".  In  the  following,  let  us  examine  what  can  be 
obtained  from  the  theory  available  in  transport  theory. 
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Fante  (1973)  has  solved  the  transport  equation  under  the  forward 
scattering  approximation,  and  Ishimaru  (1978,  ch.  13)  has  a  lucid  derivation 
and  discussion  on  it.  Here  we  only  draw  some  main  threads  for  understanding. 
Since 

A 

dI(r,Q) 

A  A 

- — -  -  2»grad  I(r,Q),  (4.1) 

dX  - 

A 

where  dX  is  the  length  of  an  elementary  segment  in  Q  direction  (Fig.  3.1),  the 
transport  equation  (3.1)  can  be  written  as 


Q*gradI(r,Q)  -  -  rieI(r,Q)+na  /  D(Q ,Q0)I(r,Q0)dQ0+w(r,Q)  .  (4.2) 

—  —  4it  — 

Because  the  scattered  energy  is  mostly  confined  within  a  small  angle  in  the 

forward  direction,  we  choose  the  z-axis  of  the  cartesian  coordinates  as  this 

direction,  and  approximate  (4.2)  through  the  following  steps. 

A  A  A  A 

Q  m  lx  +  my  +  nz  ,  (4.3) 

where  x,  y  and  z  are  the  unit  vectors  in  x,  y  and  z-axis  respectively,  and 
l,m,n,  the  corresponding  direction  cosines.  In  the  spherical  coordinate 
system  with  z-axis  as  its  polar  axis  (Fig.  4.1) 

1  ■  sinGcos#,  m  -  sinGsind,  n  -  cosG  .  (4.4) 


Because  the  angle  with  z-axis  G  is  always  small,  we  have  approximations 


tv-'  --T  -■ 


where 


_r  ■  xx-fyy+zz  ■  £  +  zz , 

*  a  „  a  A  a  * 

r  -  lx-Hny ,  Vt  -  x+^  y  . 


(4.6) 


Mote  that  £  is  not  a  unit  vector.  Because  6  is  a  small  angle,  the  magnitude 
of  _s  is  much  smaller  than  1. 

By  these  approximations  (4.2)  becomes 

31  I(s,g,s)  +  s«VtI(z,g,s) 


-tieI(zfp,s)-»—  //  D(s-s' )I(z,g,s' )ds'+W(z,g,s) 


(4.7) 


here  D(Q,Q0)  is  assumed  only  as  a  function  of  Q-Q0.  Since  most  of  the  energy 
is  confined  within  a  small  angle  with  z-axis,  the  integration  limits  for  1  and 
m  are  extended  to  l**  without  introducing  any  significant  change. 

Again  (4.7)  can  be  solved  by  the  Fourier  transform  method  (Fante  1973, 
Ishimaru  1978,  ch.  13),  the  general  solution  for  U(z,£,s)  *  0  is 

I(z,g,s)  m  I  d*  /  exp(-ik»g-is»g)I0(k,g+kz)K(z,k,g)  ,  (4.8) 


where 


Io^*3^  “  //  I0(pis)exp(ik»g+is»g)dgds 


(4.9) 


is  the  double  Fourier  transform  of  the  incident  intensty  Io(£.»£.)  at  z*°»  and 


*  B0 

K(z,k,g)  -  exp{  /  T)e[  1-r—  D(g+k(z-z' ) )  ]dz' } 


(4.10) 


1  4*  .  •  «  •  4  *  .  ■  , 


*’•  .*•  1 


D(jj)  ■  //  D(s)  exp  (is  •£)<!£  . 


(4.11) 


There  Is  no  general  explicit  expression  for  (4.8)  for  a  general 
scattering  directivity  D(s).  If  we  approximate  the  strong  forward  scattering 
pattern  by  a  Gaussian  function, 

D(s)  *  4£exp(-5a2)  (4.12) 

where  (  is  a  parameter  proportional  to  (Jl0/\)2,  and  Jt0  is  correlation  length 
of  the  random  medium,  X  is  the  wavelength,  substituting  into  (4.11)  and  (4.10) 
yields 

0(3)  -  4xexp(-^-)  ,  (4.13) 


z  2 

K(z,k, 3)  -  exp{  /  T)e{ l-B0exp(—^=-) ]dz' )  .  (4.14) 

“  o  44 

Since  most  of  the ‘energy  is  confined  within  a  small  cone  along  z>axis,  we 
consider  the  case  of  a  plane  incident  wave 

I 0(£,s>  -  I06(jL>  ,  (4.15) 

Io(k.^)  -  (2*)2I05(k).  (4.16) 


From  (4.8)  we  have 
I„ 


2 

I(z,p,s)  - -  /  dk  /  d3  exp(-ik»p-is«q)5(k)exp[-T|ez+nsz  exp(-^-)  ] .  (4. 17) 

(2it)2  '  '  45 


When  the  scattering  distance  is  large,  i.e.  T|8z>>1,  the  main  contributions  to 
the  integral  in  (4.17)  come  from  the  integrands  with  small  q's.  We  can  set 

exP(~4|-)  "  45  *  (4.18) 

Therefore 


2  f 0  p j) 


lot  rs2 

---  expl-n.z  -  — ]  . 


(4.19) 


[vv 


E(z,£>  -  j(*,£)  -  /  I(z,£,a)d£  -  I0e~\z  . 


(4.20) 


We  see  thst,  under  forwsrd  scattering  approximation,  the  energy  density  decay 
with  distance  is  only  due  to  the  absorption.  That  is  because,  in  the 
approximation,  we  neglect  the  backscattering  and  the  path  length  differences 
between  the  direct  path  and  the  multiple  scattering  paths  by  letting  cos9al. 
In  Fig.  3.11  and  Fig.  3.12  we  plot  the  energy  flux  J(r)  of  strong  forward 
scattering  vs.  that  of  the  isotropic  scattering.  If  we  consider  the 
lengthening  of  travel  paths  by  multiple  forward  scattering,  the  decay  curve 
could  be  somewhere  between  these  two  extremes. 

(4.19)  gives  the  angle  distribution  of  intensities.  The  incident  wave 
has  only  intensity  in  z-direction,  after  scattering  by  the  medium,  the 
intensities  with  different  directions  have  a  Gaussian  distribution  and  the 
width  of  the  angle  distribution  broadens  with  distance.  The  loss  due  to  the 
scattering  of  energy  to  other  directions  is  compensated  by  the  gain  of 
scattered  energy  from  other  directions.  Therefore  there  is  no  energy  loss 
except  absorption.  However,  in  order  to  calculate  the  real  energy 
attenuation,  we  have  to  take  the  backscattered  energy  into  account.  Wu 
(1982a,  see  appendix  C  )  uses  a  simple  renormalization  procedure  and  sums 
up  all  the  energy  scattered  into  the  back  halfspace  as  the  energy  loss.  This 
procedure  is  similar  to  DeUolf's  "Cumulative  Forward-Scatter 
Single-Backscatter  Approximation"  in  calculating  the  backscattering  strength 
(DeWolf  1971).  Since  the  backscattered  energy  is  much  smaller  than  the 
forward  scattered  energy,  the  second  backscattered  energy  (from  the  backward 
direction  into  the  forward  direction)  is  one  order  smaller  than  the  single 
backscattered  energy.  Therefore  the  single  backscattering  loss  with  the 
renormalization  of  the  total  forward  energy  could  be  a  reasonable 
approximation  of  the  scattering  attenuation  for  the  harmonic  wave  field. 


Fro*  the  above  analysis,  in  the  case  of  strong  forward  scattering  due  to 
large  scale  inhomogeneities,  the  shape  of  the  energy  density  decay  curve  is 
insensitive  to  the  medium  albedo  B0  and  the  separation  of  scattering 
attenuation  from  absorption  becomes  more  difficult.  However,  because  the 
scattering  loss  is  much  smaller  than  the  isotropic  case,  we  can  have  some 
constraint  on  the  possible  scattering  attenuation  from  the  strength  of 
inhomogeneities.  The  shape  of  the  seismogram  evelope  in  time  domain  can  also 
give  constraints  on  the  possible  values  of  albedo  BQ.  We  will  discuss  this 
later  in  this  paper. 


S.  SEISMIC  WAVE  SCATTERING  AND  ATTENUATION  IN  HINDU  KUSH  REGION 

In  this  section  we  will  calculate  the  energy  density  distribution  with 
travel  distances  for  the  small  earthquakes  in  the  Hindu  Kush  region.  The  data 
used  are  from  the  digital  recordings  from  two  stations  in  that  area.  Between 
11  June  and  13  July,  1977,  11  smoked  paper  recorders  and  4  digital  event 
detector  recorders  were  operated  around  the  Hindu  Kush  Mountains  of 
Northeastern  Afghanistan.  The  organization  and  operation  of  the  field  work  as 
well  as  the  seismicity  and  tectonics  of  that  region  are  described  by  Roecker 
(1981),  Chatelain  et  al.  (1980)  and  Roecker  et  al.  (1982).  Fig.  S.l  is  the 
map  view  of  the  earthquake  distribution  and  the  station  locations.  In  Fig. 

3.2,  the  events  are  divided  into  groups  with  SO  km  depth  Intervals.  The 

cn 

digital  numbered  events  were  recorded  digitally ^magnetic  tapes,  which  have 
been  used  by  Roecker  et  al.  (1982)  to  calculate  the  coda  Q  and  S  wave  Q  in 
that  region  using  Akl's  single  station  methods.  Table  5.1  lists  these  events. 
We  will  use  some  of  those  events  to  calculate  the  energy  distribution  along 
the  travel  path. 

The  digital  event  recorders  were  of  the  event  detector  type  (for  details 
see  Prothero  1976).  When  the  received  signal  exceeded  the  pre-set  level  the 
recorders  were  triggered  to  record  the  event  on  magnetic  tapes.  The  buffer  of 
the  instruments  also  allowed  us  to  record  one  second  data  proceeding  the 
triggering  signal.  Each  digital  station  had  four  seismometers,  three 
components  with  high  gain  and  a  low-gain  vertical  component.  The  natural 
period  of  the  seismometers  was  4  seconds.  The  preamplifier  had  a  gain  20  db 
or  40  db  (low  gain  or  high  gain).  The  amplifier  had  a  gain  52  db  or  58  db, 
with  a  3  pole,  low-pess,  antialiasing  filter  having  a  corner  frequency  of  32 
ha.  The  response  of  the  whole  system  is  shown  in  Fig.  5.3.  After 
amplification  the  signal  was  digitized  at  128  samples  per  second,  multiplexed 


and  than  recorded  if  the  recorder  was  triggered.  The  events  recorded  usually 
were  greater  than  magnitude  3,  with  the  exception  of  a  few  close  earthquakes, 
due  the  pre-set  trigger  level. 

Because  there  are  only  a  few  stations,  it  is  difficult  to  get  the  energy 
density-distance  relation  from  a  single  event.  We  will  use  a  single  station 
method.  The  seismograms  of  different  events  with  different  distances  from  the 
station  will  be  Fourier- transformed  to  get  the  spectral  density  of  the  energy 
density  E(r)  for  the  corresponding  source-sensor  distances.  In  order  to  have 
a  common  source  factor  for  all  the  events,  we  use  the  coda  spectral  density  of 
these  events  as  the  reference  levels.  From  observations,  it  is  generally 
acknowledged  that  the  coda  level,  at  the  travel  time  greater  than  twice  the 
S  wave  travel  time,  has  a  very  stable  relation  with  the  source  energy  and  does 
not  change  with  the  location;  of  the  events.  This  can  be  explained  by  the 
theory  of  coda  generation  in  which  the  coda  waves  are  assumed  to  be  formed  by 
the  backscattered  S  waves  from  the  heterogeneities  in  the  local  region  of  the 
lithosphere  (Aki  1969,  Aki  and  Chouet  1975).  A  received  signal  can  be 
considered  as  a  product  of  three  factors: 

received  signal  -  source  factor  x  path  factor  x  station  factor.  (5.1) 
Because  the  coda  energy  at  a  specified  time  interval  is  assumed  to  be  the  sum 
of  backscattered  wave  energy  from  the  heterogeneities  in  all  the  directions, 
therefore  the  path  factor  has  been  averaged  over  all  the  directions,  which  is 
much  more  stable  than  the  path  factor  of  the  direct  path. 

In  the  calculations,  we  took  the  reference  coda  travel  time  as  tQ  «  70 
sec.  However,  for  the  very  close  events,  some  seismograms  are  shorter  than  70 
sec,  while  for  the  distant  events,  70  sec  is  smaller  than  twice  the  S  wave 
travel  times.  We  need  to  do  extrapolations.  The  guideline  for  choosing  coda 
time  tc  is  to  have  it  greater  than  twice  the  S  travel  time  and  as  close  as 


possible  to  70  sec.  In  order  to  convert  the  coda  level  of  each  tc  to  the 
reference  level  of  tQ  ■  70  sec,  we  use  the  empirical  averaged  coda  envelope 
decay  for  each  frequency  obtained  by  Roecker  (1982)  for  this  region.  When 
t>2ts,  where  t3  is  the  S  travel  time,  the  coda  envelope  decay  can  be  fitted 
by 

1 

P(a)jt)  -  P0(o>)  —  exp(-btt),  (5.2) 

t2 

where  P(u|t)  is  the  coda  power  spectral  density  at  frequency  cu,  at  time  t, 
P0(b>)  is  a  constant,  bt  is  the  attenuation  rate  and 

bt  -  pb,  (5.3) 

where  b  is  the  attenuation  coefficient  and  p  is  S  wave  velocity.  For  the 
single  backseat tering  model,  P0(u>)  is  found  to  be  (Aki  and  Chouet  1975) 

2g(ic)S(cu) 

Po(0))  -  - ,  (5.4) 

P 

where  g(w)  is  the  backscattering  coefficient  and  S(w)  is  the  source  power. 

For  our  purpose,  it  is  not  necessary  to  specify  P0(u>),  we  need  only  use  the 
empirical  relation  (5.2).  If  we  set  t  ■  tc  as  the  reference  coda  travel  time, 
then 

1 

P((ojt0)  -  P0(u)  exp(-btt0)  .  (5.5) 

to2 

Suppose  we  measure  the  coda  power  P(u>jtc)  at  time  tc,  the  correction  for 
reducing  P(u>|tc)  to  P(u|t0)  is  then 

tc 

P(u>|t0)  »P(ujtc)  (  — )2  exp[-bt(t0-tc)  ]  .  (5.6) 

*o 

We  can  also  use  P0(io)  as  the  reference  level: 

P0(oj)  -  P(u)J  t0)  •  t02exp(btt0)  .  (5.7) 

In  Fig.  5.4,  the  solid  line  is  the  averaged  attenuation-frequency  relation 
obtained  by  Roecker,  the  dotted  line  is  the  smoothed  version  being  used  for 


calculations. 


We  choose  two  stations  PEN  and  CHS  (Fig.  5.1),  because  there  were  many 
close  events  for  both  stations  to  confine  the  energy-distance  curves.  In 
Table  5.2  and  5.3  the  events  used  for  calculations  are  listed  in  the  order  of 
distances.  The  events  were  located  using  the  arrival  times  on  smoked  paper 
records. 

To  calculate  the  spectral  density,  we  use  the  fast  Fourier  transform 
algorithm,  and  average  the  spectral  densities  over  the  specified  bandwidths. 

In  order  to  compare  with  the  previous  results  obtained  using  the  filtering 
method  by  other  authors,  we  take  the  frequencies  as  octave  and  with  bandwidths 
2/3  of  the  central  frequencies.  Table  5.4  lists  the  14  central  frequencies 
and  their  corresponding  attenuation  values.  We  use  a  32  second  window  for  the 
S  wave  Fourier  transforms.  Fig.  5.5  shows  some  examples  of  the  seismograms  at 
station  PEN  for  different  hypocenter  distances,  from  which  we  can  see  that  the 
32  second  window  will  include  most  of  the  S  wave  energy.  In  the  figure,  for 
each  event  first  gram  is  the  low  gain  vertical  (Z)  component,  the  rest  are 
high  g.  in  Z,  E-U,  and  N-S  components  respectively.  In  order  to  avoid  the 
Gibbs  phenomena  of  the  rectangular  window,  we  use  a  1  second  cosine  taper  for 
both  edges  of  the  window.  For  the  reference  coda  spectrum,  we  use  an  8  second 
Hamming  window  for  Fourier  transforms. 

Fig.  5.6  shows  the  obtained  4nr^E(r)  curves  from  the  station  PEN. 

Totally  31  events  are  used  and  the  events  are  grouped  according  to  their 
distances.  From  left  to  right,  the  curves  are  of  Z,  EW  and  NS  components.  In 
the  upper  part,  they  are  for  f  -  0.25,  0.5  and  l  kz;  in  the  middle,  f  -  1.5-8 
hz;  in  the  bottom,  f  “  12-45  hz.  Except  for  the  low  frequencies  f<l  hz ,  the 
curves  can  almost  be  fitted  with  straight  lines.  We  calculated  the  apparent 
attenuations  for  different  frequencies  for  the  EW  components  and  listed  in  the 


Table  5.5 


Because  of  the  fluctuations  of  the  measured  curves  and  the  insensitivity 
to  albedo  B0  when  Bo<0.5,  we  can  not  determine  exactly  the  values  of  B0  for 
each  frequency.  However,  we  can  get  some  constraints  on  the  B0  values  from 
the  comparison  between  the  measured  and  the  theoretical  curves. 

Since  Aki  introduced  the  single  station  method  using  S-coda  ratios  to 
measure  the  apparent  attenuations  of  short  period  body  waves  (Aki,  1980a), 
various  attenuation  mechanisms  have  been  examined  to  interpret  the 
observations,  especially  the  frequency  dependence  of  the  apparent 
attenuations.  After  discussing  different  attenuation  mechanisms,  Aki  proposed 
two  most  promising  candidates:  thermoelasticity  and  scattering  (Aki  1980a). 
However,  it  seems  only  the  scattering  mechanism  survived  in  the  literature. 
Dainty  (1981)  proposed  a  scattering  model  with  a  consta^n;  Q  medium  and 
attributed  the  observed  attenuation  as  the  sum  of  the  intrinsic  attenuations 
and  the  single  scattering  coefficient.  Assuming  an  intrinsic  Qj-2000,  he 
matched  the  observed  data  in  Kanto,  Japan  by  Aki  (1980a)  well  with  the 
theoretical  calculations.  Let  us  test  this  model  usign  our  theoretical 
calculatuions  and  the  data  in  Hindu-Kush.  Fig.  5.7  gives  the  possible  energy 
density  distribution  curves  for  different  frequencies  if  we  assume  the 
constant  Q  model  (Q^  «  2000)  and  use  the  values  of  apparent  Q  in  Kanto  region 
obtained  by  Aki.  Due  to  the  low  intrinsic  attenuation  at  low  frequencies,  the 
medium  albedo  B0  will  be  very  high,  if  we  attribute  the  observed  apparent 
attenuations  mainly  to  scattering.  However,  from  Fig.  5.7  and  Fig.  3.6,  we 
see  that  the  E^r)  curves  for  B0  1  are  o'.  arch  shape,  only  approach 
approximately  exponential  curves  when  distances  are  much  greater  than  the 
extinction  length  Le.  Compare  the  prediction  of  Fig.  5.7  with  Fig.  5.6,  they 
do  not  agree  in  general.  More  detailed  comparison  is  shown  in  Fig.  5.8  for 
the  Hindu-Kush  data.  The  apparent  attenuations  b  obtained  from  the  curves 
in  Fig.  5.6  are  listed  in  Table  5.5.  For  the  highest  frequency  f  -  45  hz, 


b  "  0.03/km.  If  we  assume  this  is  totally  due  to  the  intrinsic  absorption, 
the  equivalent  will  be  around  2500.  From  this  Q^,  we  can  obtain  the 
approximate  B0,  d0,  and  Le  for  each  frequency  (also  listed  in  Table  5.5)  based 

A 

on  the  measured  b  values.  In  Fig.  5.8  the  prediction  of  the  constant  Q  model 
with  B0  "  0.9  is  compared  with  the  measured  data  of  f  -  1.5  hz  and  2  hz. 

There  is  no  match  between  them.  Comparing  with  other  theoretical  curves  of 
different  B0's,  we  found  that  the  observations  may  be  fitted  with  curves  of 
Bo<0.5.  Since  the  fluctuation  in  the  data  and  the  insensitivity  of  the  curve 
to  B0  for  Bo<0.5,  we  can  not  determine  the  B0  value  precisely.  For  the 
frequencies  above  1.5  hz,  we  can  have  the  similar  conclusions  (see  Fig.  5.6). 
Therefore,  the  constant  model  may  not  represent  the  real  medium  in 
Hindu-Kush  region.  Because  the  apparent  attenuations  and  their  frequency 
dependence  for  some  tectonically,  active  regions  (e.g.  California,  Kanto  region 
of  Japan,  etc.)  are  alike,  we  might  expect  that  these  analyses  would  be 
applicable  to  that  region.  However,  we  need  to  apply  the  method  to  other 
regions  before  we  can  draw  any  conclusions. 

More  careful  studies  are  needed  for  the  energy  distribution  curves  of  low 
frequencies  ( f < 1  hz).  The  curves  at  these  frequencies  (Fig.  3.6)  have  the 
interesting  arch  shapes,  which  might  indicate  the  existence  of  strong 
scattering  at  these  frequencies.  However,  these  curves  are  more  fluctuating, 
which  may  be  caused  by  the  interference,  and  therefore  are  less  reliable. 
Another  consideration  is  the  influence  of  surface  waves  and  guided  waves 
(higher  mode  Rayleigh  and  Love  waves),  which  is  stronger  at  low  frequencies. 
Fig.  5.9  shows  two  examples  of  seismograms  of  events  having  distances  about 
100  km  from  the  station  (A34:  r  «  104  km,  depth  ■  4.57  km;  A08:  r  *  124  km, 
depth  16.27  km).  The  strong  low  frequency  components  following  immediately 
the  S  arrivals  are  apparent.  This  may  evidence  the  strong  multiple  scattering 


t  these  low  frequencies.  Another  positive  indicate of  multiple  scattering  is 
emerged  when  we  compare  with  the  decay  curves  of  direct  S  amplitudes.  In  Fig. 
5.10  these  decay  curves  are  shown  for  f  ■  0.25,  0.5,  and  1  hz,  the  S  wave 
power  spectra  are  calculated  by  Fourier  transform  using  a  4  second  Hamming 
window.  These  curves  are  more  regular  and  are  not  of  arch  shape.  That  is 
because  when  the  window  for  S  wave  is  very  narrow,  the  multiple-scattered 
waves,  which  have  longer  travel  times,  are  not  included.  In  5.11  we  also  plot 
the  calculated  apparent  attenuations  from  both  the  direct  S  and  the  total  S 
decay  curves  for  comparisons  (the  smoothed  coda  attenuation  curve  is  also 
plotted).  Above  1.5  hz,  the  attenuation  of  the  total  S  wave  is  smaller  than 
that  of  the  direct  S  wave.  This  may  be  due  to  the  inclusion  of  part  of  the 
scattered  energy  in  the  former  case.  However,  the  differences  between  these 
two  cases  are  small  in  general,  which  further  suggest  that,  the  scattering  is 
not  the  dominant  factor  in  the  apparent  attenuation  for  these  frequencies. 
Again  a  noticeable  different  behavior  at  low  frequencis  (f<1.5  hz)  is 
presented.  For  these  frequencies,  the  attenuations  of  direct  S  waves  are 
smaller  than  that  of  the  total  S  waves.  Note  that  the  attenuations  of  the 
total  S  waves  are  estimated  from  only  the  later  part  of  the  energy 
distribution  curves. 

If  we  take  the  energy  curves  for  f<l  hz  as  controlled  by  multiple 
scattering.  A  rough  estimation  by  comparing  with  the  theoretical  curves  (Fig. 
3.6)  can  be  made  about  the  medium  scattering  parameters.  For  the  vertical 
component,  we  have 
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In  the  case  of  the  EW  component, 


f 

Le(km) 

B0 

Ls(km) 

La(km) 

Q1(equii 

0.25 

25 

0.99 

25.3 

2500 

1122 

0.5 

26.7 

0.99 

26.9 

2670 

2396 

1 

40 

0.90 

44 

400 

718 

It  is  interesting  to  note  that,  at  0.25  and  0.5  hz,  the  estimated  are 
close  to  the  proposed  intrinsic  Q  for  the  constant  Q  model.  Although  this  may 
be  only  a  numerical  coincidence,  we  would  like  to  report  it  here  for  further 
study. 

Fig.  5.12  shows  the  4  r  E(r)  curves  obtained  from  station  CHS.  The. 


events  used  are  listed  in  Table  5.3.  The  general  conclusions  drawn  from  the 
analysis  of  the  results  of  station  PEN  hold  true  also  for  CHS. 


6.  DIFFUSION  APPROXIMATION  IN  TIME  DOMAIN,  THE  CONSTRAINT  OF  SEISMOGRAM 
ENVELOPE  ON  THE  SCATTERING  STRENGTH 

Another  approach  for  studying  the  scattering  and  attenuation  of  seismic 
waves  is  to  formulate  the  problem  of  energy  transfer  in  the  time  domain  and 
compare  the  envelopes  of  seismograms  with  the  theoretical  predictions. 

However,  from  the  author's  knowledge,  the  complete  solution  of  energy  transfer 
in  time  domain  is  not  available  at  present.  Nevertheless,  there  are 
approximate  solutions  for  the  weak  scattering  and  the  strong  scattering.  In 
the  weak  scattering  case,  when  the  propagation  distance  is  smaller  than  the 
scattering  mean  free  path,  the  single  scattering  approximation  can  be  used. 

Aki  and  Chouet  (1975)  developed  a  single  backscattering  model,  Sato  (1977) 
derived  the  formulation  for  isotropic  scattering  and  discussed  subsequently 
the  influence  of  non-isotroplc  scattering  (Sato  1982).  In  the  case  of  strong 
scattering,  when  the  scattering  coefficient  is  much  greater  than  the 
absorption  coefficient  (Bo>>0.5),  and  the  propagation  distance  is  much  greater 
than  the  scattering  mean  free  path,  the  diffusion  approximation  can  be  used  to 
approximate  the  envelope  variation  in  time  domain.  In  the  following,  we 
discuss  the  diffusion  approximation  and  seek  the  constraint  of  the  observed 
envelopes  on  the  medium  scattering  properties. 

When  the  scattering  mean  free  path  is  much  shorter  than  the  absorption 
length  in  the  medium,  the  energy  transfer  can  be  approximated  by  a  diffusion 
equation  (see  Morse  and  Feshbach,  1953,  §2.4) 

8 

—  P(r,t)  -  dV2P(r,t)  -  btP(r,t)  +  q(t),  (6.1) 

dt 

where  P(r,t)  is  the  power  at  distance  r  and  time  t;  bt  is  the  absorption  rate 

bt  -  be, 

where  b  is  the  absorption  coefficient  and  c  is  the  wave  velocity;  q(t)  is  the 
source;  and  d  is  the  diffusivity 
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(6.3) 


where  114  is  the  effective  extinction  coefficient  for  the  diffusion  process. 

In  the  case  of  isotropic  scattering  1^4  -  T)e.  For  non-isotropic  scattering,  in 
the  case  of  discrete  random  media 

’Id  “  afla  +  nom»  (6.4) 

where  n  is  the  number  density  of  the  scatterers,  aa  is  the  absorption 
cross-section  of  the  scatterers,  <Jm  is  defined  by 

A 

°m  “  /  "  cos9)dQ  ,  (6.5) 

4x 

A 

where  o^(Q)  is  the  differential  scattering  cross-section  (2.15),  9  is  the 
seatterng  angle.  om  is  called  the  "momentum  transfer  cross  section"  by  Morse 
and  Feshbach  (1953,  p.  188).  The  solution  of  (6.1)  for  a  point  Impulsive 
source  is  (Morse  and  Feshbach,  1953) 


P(r, t)  -  | 


(4xdt)3/;! 


-  expl-(r2/4dt)-btt] ,  t>0  . 


(6.6) 


Ishimaru  (1978)  formulated  the  problem  using  the  equations  for  the 
two-frequency  mutual  coherency  functnion  and  derived  the  two-frequency 
equation  of  transfer.  Under  the  diffusion  approximation,  a  solution  similar 
to  (6.6)  for  a  point  impulsive  source  was  obtained 


P(r,t)  - 


exp [ -(r  /4dt)-btt], 


(6.7) 


r4itd  t3/2 


where  d  is  the  same  as  (6.3),  but  with  T)<j  defined  by 


T) d  “ 


(6.8) 


instead  of  (6.4).  However,  since  aa  is  assumed  very  small,  there  will  be  no 
big  difference  between  (6.8)  and  (6.4)  except  for  very  strong  forward 
scattering.  In  the  following  we  will  discuss  the  case  of  0m>>0a>  therefore 


(6.8)  will  be  used,  which  can  be  written  as 

T>d  -  ns  (6.9) 

where  17,  is  the  scattering  coefficient,  and  y  is  the  mean  scattering  angle 
cosine 

1 

y  ■  —  /  ad(Q)cos0  dQ 
as  4it 

1 

-  —  /  D(Q)cos0  dQ,  (6.10) 

4it  4it 

A 

where  D(Q)  is  the  scattering  directivity  (2.11). 

Note  that,  the  quantity 

Tli  2*  * 

17 sY  -  —  /  d$  /  D(9  $)cos0  sin0  d0,  (6.11) 

o  o 

where  9,  $  consistute  the  spherical  coordinates  with  the  polar  axis  in  the 
incident  direction,  is  the  net  scattering  power  flux  in  the  incident 
direction.  This  part  of  the  scattered  power  will  join  the  incident  power  flow, 
and  does  not  contribute  to  the  diffusion  process.  In  the  case  of  isotropic 
scattering,  y  -  0,  the  net  scattering  power  flow  in  the  incident  direction  is 
zero,  so  that  ■  t)s.  In  the  case  of  strong  backscattering,  -l<y<0,  so 
■nd>Tls*  Vice  versa,  for  the  case  of  strong  forward  scattering,  0<y<l,  Tid<ris. 

From  (6.6),  we  know  there  is  a  peak  in  the  power  flow  curve,  which  is 

approximately  at  the  maximum  of  the  exponent  of  the  exponential  term,  l.e.  at 
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/3  _  B0 

2  1-y  l-B0 


(6.12) 


where 


are  the  absorption  time  and  the  scattering  mean  free  time  respectively. 
Therefore,  the  arrival  time  of  the  peak  of  the  power  flow  is  proportional  to 
the  square  root  of  the  ratio  between  the  absorption  time  and  the  scattering 

mean  free  time.  For  strong  forward  scattering,  /1-y  <  1,  the  power  peak  will 
arrive  earlier  than  the  case  of  isotropic  scattering;  in  the  case  of  strong 

backscattering,  /1-y  >  1,  the  peak  will  arrive  later. 

Note  that,  under  diffusion  approximation,  the  apparent  attenuation,  when 
t  xa»  t02 ,  is  approaching  to  the  absorption  coefficient  t)a;  while  in  the 
exact  solution  in  frequency  domain  (3.16)  it  approaches  or  T)a  +  d3r|s. 

The  multiplier  d0  or  da  varies  depending  on  B0.  Only  in  the  case  of  B0+l,  the 
apparent  attenuation  approaches  T)a. 

From  the  peak  time  we  can  derive  the  ratio  T)<j/na  after  doing  correction 
3/2 

of  t  ,  while  measuring  apparent  attenuation  will  determine  approximately  r|a. 
Therefore  the  shape  of  the  envelope  provides  all  the  parameters  of  diffusion 
scattering. 

If  we  assume  a  constant  Q  model,  from  table  5.5,  we  have  Bo>0.79  for  f<6 
hz.  Therefore,  the  diffusion  approximation  could  be  applied  to  the  wave 
energy  transfer  for  frequencies  below  6  hz.  Based  on  the  estimated  scattering 
parameters  in  Table  5.5,  we  list  in  Table  6.1  the  predicted  arrival  times  of 
peak  power  for  different  frequencies.  Except  for  the  strong  forward 
scattering  case,  the  peak  arrival  times  have  a  large  delay  up  to  several  times 
of  the  direct  travel  time.  This  contradicts  the  observations  on  earthquake 
seismograms  or  explosion-source  seismograms  on  the  earth.  The  travel  time 
fluctuations  for  local  earthquakes  are  usually  less  than  10-202  and  the  direct 


S  waves  can  be  easily  recognized  for  these  frequencies  in  general.  The 
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observed  seismograms  are  not  of  diffusion  type  in  the  frequency  range  1-10  hz. 
Fig.  6.1  shows  the  diffusion  type  envelope  curves  according  to  (6.6)  at  f  *  2 
hs  for  different  scattering  patterns  (the  envelopes  should  be  symmetric  about 
the  time  axis).  Since  we  neglected  the  t“  term  in  deriving  (6.12),  the 
peak  times  in  Fig.  6.1  are  different  from  the  predictions  in  Table  6.1. 
However,  the  envelopes  exhibit  the  typical  diffusion  characteristics.  These 
diffusion  type  envelopes  have  been  observed  on  the  moonquake  seismograms  and 
on  the  seismograms  of  model  experiments  in  laboratories.  In  Fig.  6.2,  the 
3-component  seismograms  for  two  events  on  the  moon  are  shown  (the  figures  are 
from  Latham  et  al.  1971).  The  first  event  (upper  seismograms)  is  believed  to 
be  a  meteoroid  impact,  corresponding  to  the  case  of  shallow  source;  the  second 
event  is  considered  to  be  a  deep  moonquake  (below  the  strong  scattering 
layer).  These  diffusion  type  seismograms  are  due  to  the  existence  of  the  high 
Q,  strong  scattering  layer  below  the  moon  surface  (Dainty  et  al.,  1974,  Dainty 
and  Toksoz,  1981).  Fig.  6.3  shows  the  seismograms  from  the  model  experiment 
in  the  laboratory  (Dainty  et  al.,  1974).  (a)  is  the  seismogram  with  a 

homogeneous  plate  as  the  propagation  medium;  while  (b)  shows  the  diffusion 
type  seismogram  for  the  case  when  the  plate  has  many  grooves  as  scatterers. 

In  order  to  compare  with  Fig.  6.1,  we  select  two  events  A06  (Depth  103 
km)  and  A15  (depth  118  km),  which  have  distances  around  200  km  from  station 
PEN  and  CHS.  From  Fig.  3.9  we  know  that,  the  diffuse  term  will  dominate  after 
the  travel  distance  exceeds  twice  the  extinction  distance  for  B0  -  0.9. 
Therefore  the  seismograms  for  these  two  events  should  be  of  diffusion  type,  if 
the  parameters  in  Table  6.1  are  true,  i.e.  the  constant  Q  model  is  true.  Fig. 
6.4  and  6.5  show  the  filtered  seismograms  for  these  two  events  at  different 
stations.  The  digital  filter  is  a  six  pole,  zero-phase,  Butterworth  filter, 
the  central  frequencies  are  0.375,  0.75,  1.5,  3,  6,  12,  24,  and  46  hz,  the 
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bandwidth  is  2/3  of  the  central  frequency  for  each  band  (see  Roecker,  1982). 
From  these  figures,  we  do  not  see  the  diffusion  type  seismograms  at  low 
frequencies.  In  Fig.  6.6  we  plot  the  envelope  decay  curves  for  event  A15  as 
an  example.  The  power  decay  curves  are  calculated  by  the  moving  window 
spectral  analysis  with  an  8  second  Hamming  window  and  at  a  3  second  interval. 
On  the  left  are  the  vertical  components,  right,  the  EW  components.  These 
envelope  curves  are  typical  for  the  events  in  this  region  (see  Roecker  1982). 
They  are  not  of  diffusion  type  except  for  some  very  high  frequency  components 
(f>20  hz ,  we  will  discuss  this  later).  In  fact  these  curves  fit  the  single 
isotropic  model  fairly  well.  The  energy  density  E(r,t)  of  the  isotropically 
scattered  body  waves  at  time  t  and  at  distance  r  from  a  point  source  can  be 
expressed  as  (Sato  1977) 


ncfsW0  t 

E(r,  t)  -  —  K(— )  , 


Aicr2 

where  ts  is  the  direct  wave  (here  S  wave)  travel  time, 
and  <Ja  the  scattering  cross-section  of  the  scatterers. 
factor,  and 


KU) 


l  5+1 

-  in  ( - ]  . 

5  5*1 


(6.13) 

n,  the  number  density 
WQ  is  the  source 

(6.16) 


The  time  function  K(t/ts)  is  a  pure  geometric  spreading  factor  for  the  single 
Isotropic  scattering  model,  which  is  plotted  in  Fig.  6.7  for  the  distance  of 
event  A15  to  CHS  (r-221.85  km).  Fig.  6i8  shows  the  power  decay  curves  after 
making  the  corresponding  geometric  spreading  correction,  i.e.  dividing  the 
curves  in  Fig.  6.6  by  K(t/ts).  We  can  see  that,  after  this  geometric 
correction,  the  power  decay  curves  are  fairly  linear,  which  is  of  exponential 
decay  due  to  attenuation. 


To  compare  with  Fig.  6.1,  w«  need  to  examine  the  case  of  strong  forward 
scatternlg  more  carefully.  The  curve  of  y*0.5  Is  calculated  by  assuming  the 
same  scattering  coefficient  J)a  as  the  isotropic  scattering  case.  Because  more 
energy  is  concentrated  in  the  forward  direction,  the  effective  scattering 
coefficient  for  diffusion  1)4  becomes  smaller  than  ti8  (see  (6.9)).  In  our  case 
we  estimated  ij8  from  the  apparent  attenuation  measurement  in  frequency  domain 
(section  5).  Since  we  calculated  the  power  spectral  densities  for  the  total  S 
waves,  the  net  scattering  power  flux  (6.11)  is  included,  so  that  the  forward 
scattering  power  flux  does  not  contribute  to  the  apparent  attenuation. 
Therefore,  the  estimated  scattering  coefficient  is  closer  to  t)<j  than  to  n8,  if 
we  consider  the  apparent  attenuatuion  is  mainly  due  to  the  scattering  loss. 

By  this  consideration,  the  curve  for  strong  forward  scattering  in  Fig.  6.1 
should  have  a  shape  close  to  the  isotropic  case  with  a  closer  to,  but  a 
*  t)8  greater  than  the  isotropic  case.  Secondly,  if  the  peak  of  the  power  flow 
is  near  the  direct  arrival  time,  the  more  elaborated  diffusion  formulae  should 
be  appealed  (Ishimaru  1978),  which  will  incorporate  the  direct  travel  time 
into  the  formulation.  At  any  rate,  if  the  apparent  attenuation  obtained  in 
section  5  is  taken  as  mainly  from  scattering  loss,  the  envelope  curve  should 
be  similar  to  a  diffusion  type  curve  of  isotropic  scattering. 

From  above  comparison  and  analysis,  combining  with  the  results  obtained 
in  section  5,  we  can  conclude  that,  in  the  frequency  range  1.5-20  hz,  the 
scattering  is  not  the  dominant  factor  of  the  measured  apparent  attenuation. 

In  other  words,  the  scattering  coefficients  is  smaller  than  the  absorption 
coefficients  at  these  frequencies  in  the  lithosphere  of  this  region. 

More  careful  study  is  also  needed  for  the  case  of  frequencies  higher  than 
20  hz.  From  Fig.  6.4,  6.5  we  notice  that,  at  these  high  frequencies  the 
seismograms  become  spindle-shaped  as  pointed  out  by  Tsujura  and  Aki  (see  Aki 


1980b).  These  are  of  diffusion  type.  For  some  stations,  the  P  and  S  phase 
can  no  Longer  be  clearly  separated,  which  means  also  strong  scattering  and 
conversion.  Since  the  attenuation  coefficients  are  high  at  these  frequencies 
the  scattering  coefficients  must  be  also  high.  This  strong  scattering  for 
high  frequencies  may  be  caused  by  the  near  surface  very  small  scale 
heterogeneities.  Regarding  Fig.  6.6,  6.8,  we  can  find  that  the  decay  curves 
of  m  and  n  band  (f-32  and  45  hz)  have  flat  tops,  different  from  the  other 
bands. 

The  time  domain  analysis  has  the  advantage  of  easy  comparison  with  the 
data,  because  each  seismogram  is  one  experiment,  unlike  the  energy  density 
distribution  curve  in  frequency  domain,  which  need  many  events  covering  a 
distance  range.  However,  in  order  to  perform  more  complete  analysis,  we  need 
to  develop  more  accurate  theory  and  model.  Besides,  the  shape  of  the  envelop 
is  also  sensitive  to  the  slip  direction  of  the  earthquake  source,  that  makes 
the  analysis  more  complicated.  At  any  rate,  the  combinations  of  time  domain 
and  frequency  domain  analysis  will  make  the  analysis  more  informative  and 


reliable. 


7.  SUGGESTION  FOR  FURTHER  STUDIES 

It  is  interesting  and  beneficial  to  apply  the  method  to  other  regions 
to  see  the  relative  importance  of  absorption  and  scattering  for  different 
regions.  Especially  the  comparison  between  the  results  for  the  tectonically 
stable  regions,  such  as  New  England  area  or  the  central  U.S.,  and  that  for  the 
active  regions  such  as  the  results  obtained  here  for  Hindu  Kush  or  that  for 
California,  will  give  us  deeper  understanding  about  scattering  and  attenuation 
as  well  as  more  information  about  the  tectonic  activities. 

Further  improvements  on  the  scattering  theory  and  modeling  are  also 
needed,  such  as  the  influences  of  the  radiation  pattern  of  the  source,  the 
finite  thickness  of  the  lithosphere,  the  nonisotropicity  of  the  inhomogeneities, 
etc.  Of  course,  full  treatment  of  elastic  wave  scattering  in  both  the  frequency 
domain  and  the  time  domain  are  highly  desired. 
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629 

636  27.30 

36 

21.49 

71. 

10. 14 

105.56 

43 

3.5 

77 

629 

1031 

4.37 

36 

28.76 

71 

18.43 

138.44 

44 

3.4 

77 

629 

1521 

33.18 

37 

27.5? 

72 

21.65 

221.95 

45 

4.3 

77 

629 

1540 

1.48 

34 

46.64 

70 

58 . 55 

2.C4 

46 

2.5 

77 

629 

16  6 

30.59 

36 

24. CO 

71 

9.33 

103.31 

06 

4.7 

77 

630 

220 

47.74 

36 

29.06 

70 

26.98 

219.25 

07 

4.1 

77 

630 

333 

34.03 

36 

37.07 

71 

17.14 

86.83 

47 

3-5 

77 

630 

1353 

29.92 

36 

17.64 

71 

11.71 

98.82 

48 

3-4 

77 

7  1 

348 

32.15 

34 

38.56 

70 

23.55 

16.27 

08 

4.5 

77 

7  1 

1444 

10.80 

36 

28.05 

71 

6.59 

264.51 

09 

4.9 

77 

7  1 

1627 

3.09 

36 

14.55 

70 

19.05 

1 Gc .83 

49 

3.0 

77 

7  2 

330 

48.31 

36 

34.45 

70 

39.30 

174.84 

50 

4. 1 

77 

7  2 

2028 

19.63 

35 

13.97 

69 

25.26 

8.32 

51 

2.2 

77 

7  2 

2111 

49.19 

35 

59.39 

70 

43.29 

93.93 

52 

3-2 

77 

7  3 

17  0 

4.91 

36 

56.43 

71 

2.82 

76.67 

53 

3.0 

77 

7  4 

614 

15.39 

36 

32.71 

71 

21.16 

122.09 

54 

3-5 

77 

7  4 

824 

3.03 

36 

19.90 

69 

33.71 

134.72 

55 

3.6 

77 

7  4 

1128 

47.08 

56 

26.39 

70 

12.86 

221.59 

10 

4.6 

77 

7  4 

2041 

5.58 

36 

11.73 

69 

26.72 

128.00 

56 

3.3 

77 

7  4 

21  1 

56.74 

37 

31.97 

72 

0.62 

157.19 

11 

4.7 

77 

7  5 

14  7 

12.46 

36 

29.72 

69 

47.10 

271.19 

57 

2.9 

77 

7  6 

055 

22.50 

36 

39.45 

71 

5.04 

229.10 

58 

3.9 

77 

7  6 

1328 

56.27 

37 

4.49 

71 

34.93 

96.33 

59 

3.4 

77 

7  6 

1659 

8.76 

36 

17.72 

69 

50.94 

9.46 

60 

3.0 

77 

7  7 

620 

43.34 

36 

25.26 

70 

37.26 

228.86 

61 

4.1 

77 

7  8 

130 

36.15 

36 

38.20 

71 

8.46 

214.61 

12 

4.0 

77 

7  8 

525 

26.60 

36 

41.86 

71 

12.24 

230.79 

13 

4.7 

77 

7  8 

7  2 

10.28 

36 

31.90 

71 

20.23 

94.37 

62 

3.3 

77 

7  8 

950 

6.99 

36 

41.21 

71 

12.71 

233.80 

63 

3.8 

77 

7  9 

1141 

13.94 

35 

30.46 

68 

52.81 

38.00 

64 

2.7 

77 

7  9 

1211 

40.56 

37 

36.71 

71 

45.78 

129.68 

65 

3.7 

77 

7  9 

1616 

39.67 

36 

28.06 

71 

12.71 

143.71 

66 

3.3 

77 

710 

028 

19.08 

35 

41.93 

68 

38.41 

3.74 

67 

3.6 

77 

710 

1347 

18.52 

35 

6.83 

69 

21.09 

17.57 

69 

2.3 

77 

710 

1612 

22.07 

35 

31 .28 

69 

13.05 

17.04 

70 

3.6 

77 

711 

11  2 

56.61 

36 

26.37 

71 

20.63 

104.40 

14 

3-9 

77 

711 

1224 

5.07 

36 

45.50 

71 

28.64 

188.43 

72 

3.6 

77 

711 

1651 

7.33 

36 

28.89 

71 

9.82 

118.43 

15 

4.2 

77 

712 

11  7 

59.16 

36 

32.53 

70 

58.27 

192.49 

73 

4.3 

77 

712 

1518 

28.32 

36 

12.21 

69 

15.96 

5.57 

74 

3.8 

77 

712 

1718 

2.32 

36 

16.59 

70 

40.23 

103.88 

75 

3.7 

Table  5.2 


Events  deed  In  the  Calculations  for  PEN  in  the  Order  of  Distances  (31  events) 


•••  Reference  Record 

Point  No.  Event  Distance  Depth  Magnitude  P  travel  S  travel  coda  time  length 
in  curve  No.  (km)  (km)  time  (sec.)  time  (sec.)  tc  (sec.)  (sec.) 


1 

-A42 

11.12 

8.68 

2.1 

7.29 

11.83 

44 

^A51 

8.32 

2.2 

7.54 

12.26 

31.5 

34 

2 

5- A69 

21.52 

17.57 

2.3 

7.54 

12.52 

31.5 

40 

^A18 

22.03 

1.49 

2.4 

9.55 

15.62 

20.5 

18 

3 

^*A70 

37.76 

17.04 

3.6 

10.10 

17.03 

34.1 

38 

>U1 

39.15 

2.5 

3.1 

11.67 

19.42 

41.7 

40 

4 

A64 

63.09 

38.0 

2.7 

11.42 

19.79 

40.4 

38 

^A23 

65.15 

1.87 

2.5 

15.04 

25.38 

42.0 

36 

5 

A67 

80.75 

3.74 

3.6 

16.86 

28.66 

60.7 

60 

6 

A20 

100.16 

77.98 

3.7 

15.94 

27.87 

54 

\  A34 

104.68 

4.57 

3.9 

19.81 

33.91 

73.8 

82 

\A35 

107.97 

9.65 

3.3 

19.67 

33.76 

68 

\74 

108.33 

5.57 

3.8 

20.16 

34.56 

80 

7 

A08 

124.88 

16.27 

4.5 

21.06 

36.39 

158 

.  A68 

175-74 

24.21 

42.61 

60 

8 

A04 

178.52 

106.25 

3.5 

24.78 

43.42 

78.7 

70 

>A39 

179.28 

100.0 

3,6 

24.9 

43.56 

64 

9 

A16 

196.42 

102.36 

3.9 

27.16 

47.47 

84.0 

70 

A75 

197.95 

103.88 

3.7 

27.36 

47.83 

64 

<*  A06 

234.53 

103.31 

4.7 

31.95 

55.72 

105.9 

134 

^415 

247.35 

118.48 

4.2 

33.41 

58.51 

107.5 

158 

11 

A03 

259.49 

136.57 

34.71 

61.11 

62 

^  A50 

295.83 

174.84 

4.1 

34.46 

61.07 

90 

.  A10 

271.40 

221.59 

4.6 

35.57 

63.1 

109.4 

258 

12 

/  AO  5* 

278.65 

218.87 

4.6 

36.52 

64.75 

234 

\A07 

279.0 

219.25 

4.1 

36.56 

64.82 

100.5 

76 

NA73 

283.65 

192.49 

4.3 

37.22 

65.96 

90 

13 

A02 

310.97 

232.49 

4.1 

40.23 

71.32 

102.1 

78 

14 

<*A13! 

329.59 

230.79 

4.7 

42.47 

75.29 

246 

NA09* 

339.32 

264.91 

4.9 

43.36 

77.06 

302 

15 

A25 

472.39 

340.70 

5.3 

58.59 

104.48 

132.2 

96 

*high  gain  records  were 

clipped 

,  only  low 

gain  vertical  component 

has  been  used. 


t: 

Point  No. 

Event 

Distance 

Depth 

Magnitude  P  travel  S  travel  Record  Length 

1 

i 

in  curve 

No. 

(kn) 

(kn) 

time  (sec.)  time  (sec.)  (sec.) 

62 


Table  5.4  The  Central  frequencies  used  end  the  corresponding  attenuation 
▼nines  of  coda  wares 


Band 

no. 

Central 

frequency 

Code  Qc 
(observed) 

Coda  Qe 
(smoothed) 

Coda  bt 
(observed) 

Coda  bt 
(smoothed) 

Coda  b 
(smoothed) 

a 

0.25 

24.0 

24.0 

6.5  x  10-2 

6.5  x  10-2 

1.86  x  10-2 

b 

0.5 

47.9 

44.2 

6.6 

7.10 

2.03 

c 

1 

83.2 

81.0 

7.6 

7.76 

2.22 

d 

1.5 

89.1 

115.4 

10.6 

8.17 

2.33 

e 

2 

107.2 

148.3 

11.7 

8.47 

2.42 

f 

3 

125.9 

211.3 

15.0 

8.92 

2.55 

8 

4 

190.5 

271.6 

13.2 

9.26 

2.64 

h 

6 

281.8 

386.8 

13.4 

9.75 

2.78 

i 

8 

446.7 

497.2 

11.3 

10.11 

2.89 

J 

12 

707.9 

708.2 

10.7 

10.65 

3.04 

k 

16 

933.3 

910.2 

10.8 

11.04 

3.16 

1 

24 

1174.9 

1296.6 

12.8 

11.63 

3.32 

■ 

32 

1698.2 

1666.6 

11.8 

12.06 

3.45 

n 

45 

2238.7 

2244.0 

12.6 

12.60 

3.60 

63. 


Table  3.5  Apparent  attenuation*  for  the  EV  coaponents  of  station  PEN 
and  the  estimated  values  of  seismic  albedo  B0's,  if  we 
assume  a  constant  Q  (>2500)  medium 


f 

-  w/p  q_1 

<Q  -  2500) 

La(km)  (for  Q-2500) 

mm 

A 

Tl*/b 

B 

Le  (km) 

0.5 

0.036  x  10-2 

2778 

1.00  x  10-2 

0.036 

-0.96 

0.5 

50 

1 

0.072 

1389 

1.38 

0.052 

0.95 

0.5 

36 

2 

0.144 

694 

1.50 

0.096 

0.90 

0.6 

40 

3 

0.215 

465 

1.60 

0.134 

0.87 

0.6 

38 

6 

0.431 

232 

2.03 

0.212 

0.79 

0.7 

34 

2 

0.862 

116 

2.50 

0.345 

0.66 

0.8 

31 

4 

1.72 

58 

2.73 

0.63 

0.37 

0.95 

35 

5 

3.23 

31 

3.00 

_l 

1.08 

0 

1 

33 

64. 


Table  6.1  The  Predicted  arrival  time  of  the  peak  power  by  the  dif fusion 
approximation  based  on  the  assumed  constant  Q(*2500)  model  and 
the  estimated  parameters  in  Table  5.5. 


f 

<hs) 

absorption 

time 

?4  (sec.) 

mean  free 
time 

ts  (see.) 

d(l-7) 

(km2/sec) 

Albedo 

Bo 

1 

Arrival  time  of  peak  power 

tjn/  tjj 

y  -  -0.5 

D 

0.5 

793.7 

14.9 

60.8 

0.96 

1 

50 

5.20 

4.24 

3.0 

1 

396.8 

10.8 

44el 

0.95 

35 

4.62 

3.77 

2.67 

2 

198.4 

12.7 

51.9 

0.90 

40 

3.18 

2.60 

1.84 

3 

132.9 

12.5 

51.0 

0.87 

38 

( 

2.74 

2.24 

1.58 

6 

66.3 

12.3 

50.2 

0.79 

34 

2.06 

1.68 

1.19 

65. 


Figure  Captions 

3.1  The  derivation  of  the  transfer  equation  for  the  specific  intensity 

A 

I(r,Q). 

3.2  The  diffuse  multipliers  dQ  and  ds  as  functions  of  BQ  (the  medium  seismic 
albedo). 

3.3  The  behavior  of  the  integrand  of  the  integral  for  the  coherent  term. 

3.4  The  normalized  energy  density  distribution  curves  4itr2E(r),  where  r  is 
the  propagation  distance  from  the  point  source.  At  the  top  are  the 
curves  of  the  diffuse  term,  at  the  bottom  are  that  of  the  coherent  term; 
in  the  middle  are  the  curves  of  the  sum  of  the  two  term.  Here  De  is  the 
numerical  extinction  distance,  Lc  -  l/rie  is  the  extinction  length  of  the 
medium,  T)a  ■  t;s  +  r|a  is  the  extinction  coefficient,  where  t;s  and  T)a  are 
the  scattering  coefficient  and  the  absorption  coefficient  respectively. 

B0  ■  h j/ (fls+na)  ia  the  medium  seismic  albedo. 

3.5  The  schematic  diagram  of  a  possible  multiple  scattering  path  compared 
with  the  direct  path.  The  hatched  shell  of  unit  thickness  will  receive 
the  energy  4xr2E(r). 

3.6  The  normalized  energy  distribution  curves  4xr2E(4)  in  the 
semi-logarithmic  scale. 

3.7  The  energy  distribution  curves  with  the  numerical  absorption  distance  Da 

A 

■  r/La,  where  La  ■  rja  is  the  absorption  length  of  the  medium,  b  is 
the  apparent  attenuation  coefficient  obtained  from  the  slope  of  the 
curve.  B0  is  the  medium  albedo. 

3.8  The  energy  distribution  curves  with  the  numerical  scattering  distance 

D*  ■  4/Ls,  where  La  ■  l/t|s  is  the  scattering  length  of  the  medium.  B0  is 


the  medium  albedo. 


3.9  The  relative  strengths  of  the  diffuse  term  E<j  and  the  coherent  term  Ec 
at  different  extinction  distances  De  ■  r/Le  for  different  medium  albedo 
B0,  where  La  -  1/t|€  is  the  extinction  length  of  the  medium. 

3.10  Seme  as  3.9,  at  different  scattering  distances  Ds  ■  r/La,  where  La  “  l/rja 
is  the  medium  scattering  length. 

3.11  The  normalized  radial  energy  flux  density  4rcr2Jr(r)  for  the  isotropic 
scatternig  case  and  the  strong  forward  scattering  case. 

3.12  Seme  as  3.11.  The  distance  is  the  numerical  absorption  distance  Ds  - 
r/La,  where  La  ■  l/ria  is  the  absorption  length  of  the  medium. 

4.1  The  derivation  for  the  case  of  strong  forward  scattering  approximation, 
z  is  along  the  forward  direction.  _r  is  the  position  vector,  £  is  the 

A 

position  vector  in  the  transverse  plan;  Q  is  the  unit  vector  in  the 

A 

scattering  direction,  and  £  is  projection  of  Q  in  the  transverse  plan. 

3.1  Hep  view  of  se'ismicity  in  the  Hindu  Kush  as  determined  by  Chatelain  et 
al.  (1980).  The  digital  stations  are  indicated  by  open  stars,  and  the 
smoked  paper  stations  by  solid  diamonds. 

5.2  Map  view  of  all  the  Hindu  Kush  seismicity  on  smoked  paper  stations, 
divided  into  50  km  depth  intervels.  Locations  of  events  recorded  on  the 
digital  recorders  are  denoted  by  numbers  used  in  Table  5.1  (from  Roecker 
1982). 

5.3  The  overall  response  of  the  digital  recorder  (from  Roecker  1981). 

5.4  The  averaged  coda  attenuation  rate  bt  ■  (3b,  where  (3  is  the  shear  wave 
velocity,  b  is  the  attenuation  coefficient.  The  solid  line  is  obtained 
by  Roecker  (1982)  for  the  shallow  events,  and  the  dotted  line  is  the 
smoothed  curve  used  in  this  paper. 


5.5  Sob*  seismograms  for  different  hypoeenter  distances  at  station  PEN.  A42 
r  ■  11.12  km,  A67,  r  *  8.75  km;  A16,  r  -  196.42  km;  A02,  r  ■  310.97  km. 

5.6  Energy  distribution  curves  4«r2E(r)  obtained  from  the  data  at  station 
PEN.  From  left  to  right:  Z,  EW  and  NS  components.  From  top  to  bottom: 
f  -  0.25-1  Hz,  f  -  1.5-8  Hz,  and  f  -  12-45  Hz. 

5.7  The  predicted  4xr2E(r)  curves  by  the  constant  Q  (>2000)  model  for 
different  frequencies,  if  the  measured  apparent  attenuation  in  Kanto, 
Japan  by  Aki  (1980a)  is  assumed  as  the  sum  of  the  absorption  coefficient 
and  the  scattering  coefficient  (Dainty  1981). 

5.8  The  comparison  between  the  observed  4icr2E(r)  for  f  ■  1.5  and  2  hz  at 
station  PEN  in  Hindu  Kush  and  the  theoretical  predictions  for  different 
B0's.  The  curve  of  BQ  >  0.9  is  the  prediction  from  the  constant  Q 
(>2500)  model,  which  does  not  match  with  the  observation. 

5.9  Examples  of  seismograms  at  station  PEN  (A34:  r  >  104  km,  depth  *  4.57  km 
A08:  r  ■  125  km,  depth  •  16.27  km),  which  show  strong  low  frequency 
coaiponents  immediately  after  the  direct  S). 

5.10  The  energy  density  curves  4xr2E(r)  for  direct  S  waves  at  f  >  0.25,  0.5 
and  1  Hz  for  station  PEN.  The  curves  are  calculated  using  a  4  sec 
Hamming  window  for  the  direct  S  arrivals.  Compare  to  fig.  5.6.  No  arch 
shape  appears  here. 

5.11  Apparent  attenuations  derived  from  the  slopes  of  the  energy  density 
curves  (Fig.  5.6)  for  station  PEN,  together  with  the  average  coda 
attenuations  and  the  direct  S  attenuations. 

0:  for  EW  component,  total  S 

X:  for  Z  component,  total  S 

A:  Z  component,  direct  S  (4  sec.  window). 


Note:  For  f<lHz,  the  apparent  attenuations  are  calculated  by  using  only 
the  last  part  of  the  curves  (Fig.  5.6). 

5.12  The  energy  density  curves  4xr2E(r)  for  station  CHS.  From  left  to  right: 
Z  component  and  EW  component.  From  top  to  bottom:  f  ■  0.25-1  Hz,  f  ■ 
1.5-8  Hz,  and  f  -  12-45  Hz. 

6.1  The  seismogram  envelopes  of  S  waves  predicted  by  the  diffusion 

approximation  for  the  case  of  f  •  2  Hz,  B0  ■  0.9  (Qj  »  2500).  y  is  the 

mean  scattering  angle  cosine  defined  by  (6.10). 

6.2  The  seismograms  of  moonquakes.  The  event  on  the  top  is  supposed  to  be  a 
meteoroid  impact;  the  bottom  event  is  believed  to  be  a  deep  moonquake 
(from  Latham  et  al. ,  1971). 

6.3  The  seismograms  from  the  model  experiment  in  laboratory  (Dainty  et  al., 
1974). 

a)  The  seismogram  with  the  homogeneous  plate. 

b)  The  seismogram  when  the  plate  has  many  grooves  as  scatterers. 

6.4  fhe  band-pass  filtered  seismograms  of  A06  (r  11  235  km,  depth  ■  103  km) 

for  the  stations  CHS,  FRA,  JOR  and  PEN.  From  top  to  bottom:  f  -  0.375, 

0.75,  1.5,  3,  6,  12,  24,  46  Hz. 

6.5  The  band-pass  filtered  seismograms  of  A15  (r  *  247  km,  depth  -  118  km) 
for  the  stations  CHS,  PEN  and  JOR.  From  top  to  bottom:  f  *  0.375,  0.75, 
1.5,  3,  6,  12,  24,  46  Hz. 

6.6  The  envelope  decay  curves  of  A15  (r  -  247  km,  depth  -  118  km)  for 
station  PEN.  From  left  to  right:  Z  component  and  EW  component.  From 
top  to  bottom:  f  ■  0.25,  0.5  Hz;  f  ■  1-8  Hz;  f  *  12-45  Hz. 


The  theoretical  envelope  decay  curve  for  the  single  isotropic  scattering 
in  a  lossless  medium  according  to  Sato  (1977).  The  envelope  decay  is  a 
pure  geometric  spreading  effect.  The  distance  between  the  source  and 
sensor  is  taken  as  247  km  (as  the  case  of  A15  to  PEN). 

The  coda  decay  curves  at  station  PEN  for  A15  after  the  geometric 
correction.  The  corrections  were  done  by  taking  the  ratios  of  the 
curves  in  Fig.  6.6  and  that  in  Fig.  6.7.  Note  that,  the  curves  for  f  ■ 
1-20  Hz  can  be  approximately  fitted  with  straight  lines,  which  means 
that,  the  scattering  at  this  frequency  range  can  be  approximated  by  the 
single  scattering  theory. 
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3.2  The  diffuse  multipliers  d0  and  ds  as  functions  of  B0  (the  medium  seismic 
albedo) . 


3.4  The  normalized  energy  density  distribution  curves  4itr^E(r),  where  r  Is 
the  propagation  distance  from  the  point  source.  At  the  top  are  the 
curves  of  the  diffuse  term,  at  the  bottom  are  that  of  the  coherent  term; 
In  the  middle  are  the  curves  of  the  sum  of  the  two  terra.  Here  De  Is  the 
numerical  extinction  distance,  Le  *  1/^le  is  the  extinction  length  of  the 
medium,  rje  -  r;,  +  t|a  is  the  extinction  coefficient. 


Da  =  r/La 


7  The  energy  distribution  curves  with  the  numerical  absorption  distance  Da 

A 

■  r/La,  where  La  ■  t)a  is  the  absorption  length  of  the  medium,  b  is  | 

the  apparent  attenuation  coefficient  obtained  from  the  slope  of  the 


curve.  Ba  is  the  medium  albedo. 
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3.8  The  energy  distribution  curves  with  the  numerical  scattering  distance 
Ds  ■  4/Ls,  where  Ls  ■  l/rjs  is  the  scattering  length  of  the  medium.  B0 


the  medium  albedo. 
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12  Same  as  3.11.  The  distance  is  the  numerical  absorption  distance  Ds 
r/Lft,  where  La  "  i/ha  Is  the  absorption  length  of  the  medium. 


4.1  The  derivation  for  the  case  of  strong  forward  scattering  approximation, 
z  is  along  the  forward  direction.  _r  is  the  position  vector,  £  is  the 

A 

position  vector  in  the  transverse  plan;  Q  is  the  unit  vector  in  the 

A 

scattering  direction,  and  s  is  projection  of  Q  in  the  transverse  plan. 
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5.1  Map  view  of  seismicity  in  the  Hindu  Kush  as  determined  by  Chatelain  et 
al.  (1980).  The  digital  stations  are  indicated  by  open  stars,  and  the 
smoked  paper  stations  by  solid  diamonds. 

5.2  Map  view  of  all  the  Hindu  Kush  seismicity  on  smoked  paper  stations, 
divided  into  50  km  depth  intervals.  Locations  of  events  recorded  on  the 
digital  recorders  are  denoted  by  numbers  used  in  Table  5.1  (from  Roecker 
1982).  — » 


Frequency  depedence  of  the  attenuation  bt 


*•4  The  averaged  coda  attenuation  rate  bj  “  Pb,  where  p  is  the  shear  wave 
velocity,  b  is  the  attenuation  coefficient.  The  solid  line  is  obtained 
by  Roecker  (1982)  for  the  shallow  events,  and  the  dotted  line  is  the 
smoothed  curve  used  in  this  paper. 


5.5  Some  seismograms  for  different  hypocenter  distances  at  station  PEN.  A i 
r  -  11.12  km,  A67,  r  «  8.75  km;  A16,  r  -  196.42  km;  A02,  r  *  310.97  km 
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5.7  The  predicted  4itr2E(r)  curves  by  the  constant  Q  (**2000)  model  for 

different  frequencies,  if  the  measured  apparent  attenuation  in  Kanto, 
Japan  by  Aki  (1980a)  is  assumed  as  the  sum  of  the  absorption  coefficient 
and  the  scattering  coefficient  (Dainty  1981). 
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A08:  r  *  125  km,  depth  ■  16.27  km),  which  show  strong  low  frequency 
components  immediately  after  the  direct  S). 


Coda  b  (smoothed) 


f  o  EW  component!  . 

S  attenuation  J  x  Z  component  J  sec  w,n  ow 
[a  Z  component  4  sec.  window 
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3.11  Apparent  attenuations  derived  from  the  slopes  of  the  energy  density 
curves  (Fig.  3.6)  for  station  PEN,  together  with  the  average  coda 
attenuations  and  the  direct  S  attenuations. 

0:  for  EW  component,  total  S 

X:  for  Z  component,  total  S 

A:  Z  component,  direct  S  (4  sec.  window). 


Note:  For  f<lHz,  the  apparent  attenuations  are  calculated  by  using  only 
the  last  part  of  the  curves  (Fig.  5.6). 
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6.1  The  seismogram  envelopes  of  S  wmi  prtdlctad  by  the  diffusion 
approximation  for  tha  ease  of  f  ■  2  Hs,  B0  -  0.9  (Q*  -  2500). 
moan  seattaring  angle  cosine  defined  by  (6.10). 
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Compressed  time-vctle  records  of  two  of  the  lunar  seismic  events  believed  to  be  of  natural 
origin  recorded  at  station  12.  Z  n  the  vertical  component  seismometer:  X  and  )'  are  the  horizontal 
component  seismometers.  The  moonquake.  event  of  13:09  hr..  May  23.  I9?0,  originated  within  the 
tone  of  greatest  activity  (At  to  net.  The  H-phase  is  prominent  on  the  seismograms  from  the  horizontal 
component  seismometers  for  category  rtt  events.  This  phase  is  tentatively  identified  as  the  direct 
shear  wave  arrival.  The  event  of  8:09  hr.,  April  8,  !9'0,  is  believed  to  he  a  meteoroid  impact 

(category  C  cventi. 


sal  sinograms  of  aoonquakas.  The  event  on  the  top  Is  supposed  to  be  a 
aeteorold  impact;  the  bottom  event  is  believed  to  be  a  deep  moonquake 
(from  Latham  et  al. ,  1971)* 
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6.3  The  (litwgKU  from  the  model  experiment  in  laboratory  (Dainty  et  al. , 
1974). 

a)  The  seismogram  with  the  homogeneous  plate. 

b)  The  seismogram  when  the  plate  has  many  grooves  as  seatterers. 


6.4  The  hand-pass  filtered  seismograms  of  A06  (r  -  235  km,  depth  -  103  km) 
for  the  stations  CHS,  FRA,  JOR  and  PEN.  From  top  to  bottom:  f  -  0.375, 
0.75,  1.5,  3,  6,  12,  24,  46  Hs. 


6.5  The  band-pass  filtered  seismograms  of  A15  (r  -  247  km,  depth  “  118  km) 

for  the  stations  CHS,  PEN  and  JOR.  From  top  to  bottom:  f  -  0.375,  0.75, 
1.5,  3,  6,  12,  24,  46  Hs. 
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The  theoretical  envelope  decay  curve  for  the  single  isotropic  scattering 
in  a  lossless  medium  according  to  Sato  (1977).  The  envelope  decay  is  a 
pure  geometric  spreading  effect.  The  distance  between  the  source  and 
sensor  is  taken  as  247  km  (as  the  case  of  A15  to  PEN). 
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